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Joel  H.  Ferziger 
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Stanford  University 
Stanford  CA 

August  1998 


1  Introduction 

Environmental  and  geophysical  flows  occur  on  much  larger  scales  than  flows  in  industrial  devices; 
one  consequence  is  that  the  Earth’s  rotation  often  plays  a  significant  role  in  them.  Furthermore, 
the  horizontal  scales  are  generally  much  larger  than  the  vertical  ones  so  that  the  motions  on 
the  largest  scales  are  essentially  two  dimensional.  Finally,  stratification  is  often  important  in 
these  flows,  sometimes  so  much  so  that,  despite  the  large  scales  and  high  Reynolds  numbers, 
the  flow  state  is  essentially  laminar.  All  of  these  effects  make  simulation  of  environmental  flows 
very  challenging  and  different  from  industrial  flows. 

One  of  the  most  difficult  issues  in  simulation  of  these  flows  is  the  need  to  deal  with  the 
real  world  in  real  time.  In  weather  or  ocean  prediction  the  goal  is  to  predict  the  actual  state  of 
the  system  at  some  time  not  too  far  in  the  future  rather  than  the  statistical  quantities  that  are 
typically  needed  in  engineering.  This  requires  simulation  of  a  single  realization  rather  than  an 
ensemble  average  flow,  providing  a  severe  constraint  on  the  approaches  that  can  be  used.  There 
is  a  limit  on  how  far  into  the  future  the  state  can  be  predicted — the  well-known  ‘butterfly’ 
problem.  (Exceptions  are  climate  studies  and  some  types  of  pollution  studies.) 

The  field  is  vast  and  no  one  paper  or  even  a  book  could  cover  all  of  it.  In  this  paper,  we 
shall  concentrate  on  some  issues  in  small  scale  meteorology  and  oceanography.  Although  they 
deal  with  the  smallest  scales  of  interest  in  these  fields,  the  dimensions  are  very  large  compared 
to  those  found  in  industrial  applications.  We  may  be  dealing  with  the  atmospheric  boundary 
layer  (typically  1-2  km  deep),  the  oceanic  mixed  layer  (typically  100-300  m  deep),  an  estuary,  or 
the  atmosphere  in  an  urban  area  and  its  surroundings;  these  have  horizontal  dimensions  of  tens 
of  kilometers.  On  these  scales,  three  dimensionality  is  important  but  Coriolis  and  stratification 
effects  may  also  play  major  roles. 

The  difficulties  of  simulating  these  flows  are  considerable  but  the  value  of  accurate  predic¬ 
tions  is  enormous;  the  effects  of  weather  and  the  state  of  the  ocean  are  constantly  in  the  news. 
We  shall  show,  that  even  in  some  of  the  simpler  flows,  there  remain  issues  to  be  resolved  that 
differ  from  those  encountered  in  industrial  flows.  Then  we  shall  give  some  examples  of  some 
recent  simulations  in  this  area. 
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2  Environmental  vs.  Industrial  Flows 


With  the  exception  of  some  very  simple  but  important  building  block  flows,  direct  numerical 
simulation  (DNS)  of  turbulent  flows  is  beyond  the  capability  of  present  computers  and  will 
remain  so  for  the  indefinite  future.  For  industrial  flows,  this  leaves  one  with  the  choice  of 
doing  large  eddy  simulations  (LES)  or  computations  based  on  the  Reynolds  averaged  Navier- 
Stokes  (RANS)  equations  which  employ  well-known  one  point  turbulence  models.  The  trade-off 
is  a  simple  one — the  low  cost  of  RANS  vs.  the  improved  accuracy  and  detail  of  LES.  In  the 
environmental  area,  the  RANS  concept  might  be  applicable  in  some  applications  such  as  climate 
the  prediction  of  the  average  state  of  an  estuary  but  weather  or  oceanic  state  predictions  must 
produce  an  accurate  estimate  of  the  actual  future  state.  For  that  purpose,  there  is  no  choice 
but  to  use  LES.  Indeed,  LES  originated  in  meteorology  as  a  response  to  this  need. 

As  already  noted,  a  typical  LES  of  the  atmospheric  boundary  layer  or  the  oceanic  mixed 
layer  simulates  a  domain  whose  horizontal  dimensions  are  hundreds  of  meters  to  tens  of  kilome¬ 
ters.  The  surfaces  might  be  idealized  as  smooth  but  they  are  almost  always  rough  on  the  scale  of 
a  typical  near-surface  turbulent  eddy.  As  these  flows  do  not  normally  separate  (separation  does 
occur  in  the  presence  of  topography)  and  because  it  is  impossible  to  represent  the  surface  (let 
alone  such  features  as  trees)  exactly,  it  is  usually  represented  by  means  of  an  artificial  boundary 
condition.  Using  such  a  condition  also  eliminates  the  need  for  a  fine  grid  near  the  surface  which 
increases  the  cost  of  engineering  simulations.  Due  to  the  size  of  the  domain  and  the  fact  that 
only  about  a  hundred  grid  points  can  be  used  in  each  coordinate  direction,  the  grid  size  or  filter 
width  must  be  tens  of  meters,  which  is  much  larger  than  the  Kolmogoroff  scale  (typically  a 
few  centimeters).  This  would  be  no  problem  for  homogeneous  turbulence  with  a  filter  scale  in 
the  inertial  subrange  but  significant  dynamics  occur  at  scales  smaller  than  the  grid  size.  For 
example,  entrainment  takes  place  in  small  zones  at  the  edge  of  a  plume  or  in  the  inversion  at  the 
top  of  the  atmospheric  boundary  layer,  a  phenomenon  that  is  difficult  to  simulate  and  remains 
a  subject  of  research. 

Consequently,  there  may  be  a  dynamic  length  scale  smaller  than  the  subgrid  scale  cutoff. 
This  makes  .it  necessary  to  use  a  turbulence  model  in  which  the  length  scale  is  not  pre-selected 
(as  in  the  Smagorinsky  model)  and  subgrid  scale  modeling  may  then  resemble  RANS  modeling 
more  than  it  does  the  kind  of  subgrid  scale  modeling  used  for  laboratory  scale  flows.  It  also 
follows  that  the  similarity  required  for  the  dynamic  procedure  to  work  may  not  be  available  in 
these  flows. 

At  the  larger  scales  of  meteorology  and  oceanography,  the  flows  become  almost  two  di¬ 
mensional  and  the  Coriolis  force,  which  plays  a  role  at  the  small  scales,  becomes  dominant. 
It  goes  without  saying  that  DNS  is  totally  impossible  at  these  scales  and  the  meaning  of  LES 
needs  to  be  reconsidered.  Solution  of  the  Navier-Stokes  equations  is  impossible  and  simplifica¬ 
tions  of  them  are  often  used;  it  is  not  unusual  to  solve  different  sets  of  equations  at  each  scale. 
The  averaging  scales  are  now  hundreds  of  kilometers  and  almost  all  of  what  an  engineer  would 
consider  turbulence  lies  in  the  subgrid  scale.  Despite  the  dominance  of  two  dimensionality,  the 
weak  motions  in  the  vertical  direction  play  an  essential  role  and  cannot  be  ignored  entirely.  At 
these  scales  the  parameterization  of  the  subgrid  scales  is  required  to  represent  a  very  wide  range 
of  phenomena  and  the  difficulty  of  constructing  them  is  extreme.  A  goal  of  smaller  scale  simu¬ 
lations  is  to  provide  these  parameterizations;  this  task  is  made  very  difficult  by  the  large  range 
of  parameters  and  conditions  needed  to  specify  the  state  of  the  smaller  scale  flow.  Furthermore, 
the  similarity  on  which  the  dynamic  subgrid  scale  model  is  based  does  not  exist  at  these  scales 
because  two  dimensional  turbulence  is  qualitatively  different  in  character  from  three  dimensional 
turbulence. 
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We  shall  describe  some  simulations  relevant  to  environmental  fluid  mechanics,  both  DNS 
and  LES,  which  have  been  done  recently.  Specifically,  we  look  at  turbulence  subject  to  stratifi¬ 
cation  and  shear  and  the  problem  of  upwelling  in  the  coastal  ocean. 

3  Simulations 

3.1  Stratified  Sheared  Homogeneous  Turbulence 

Homogeneous  sheared  turbulence  has  been  studied  from  the  earliest  days  of  turbulence  simu¬ 
lation;  shear,  of  course,  tends  to  amplify  the  intensity  of  turbulence.  The  addition  of  stable 
stratification  introduces  a  force  that  tends  to  reduce  the  intensity  of  the  turbulence  by  al¬ 
lowing  conversion  of  its  kinetic  energy  into  potential  energy.  The  interplay  of  the  two  forces 
is  what  makes  the  study  of  this  type  of  flow  interesting.  The  parameter  traditionally  used 
to  characterize  the  relative  importance  of  the  two  forces  is  the  gradient  Richardson  number 
Ri  =  ( g/p)(dp/dz)/(dU/dz )2.  Recent  studies  of  inhomogeneous  flows  show  that  this  parameter 
does  not  correlate  the  data  well;  for  example,  it  varies  by  more  than  an  order  of  magnitude  with 
distance  from  the  wall  in  stratified  turbulent  channel  flow.  A  better  parameter  is  the  turbulent 
Froude  number  Fr  =  q/NL  where  N  is  the  Brunt-Vaisala  frequency  N 2  =  (gf  p)(dp/dz).  It  can 
be  regarded  as  the  inverse  square  root  of  a  Richardson  number  but  it  is  based  on  turbulence 
rather  than  mean  flow  properties.  The  data  can  be  correlated  quite  nicely  with  this  variable. 
Especially  important  for  correlation  purposes  is  the  mixing  efficiency  which  is  the  roughly  the 
fraction  of  the  energy  extracted  from  the  mean  flow  that  goes  into  the  production  of  potential 
energy.  It  correlates  very  well  as  a  function  of  Froude  number. 

There  is  also  an  important  effect  of  the  means  of  creation  of  the  turbulence  on  the  prop¬ 
erties  of  the  resulting  turbulence.  When  shear  is  the  primary  forcing,  the  component  of  the 
turbulent  fluctuations  in  the  mean  flow  direction  (which  is  usually  horizontal  in  stratified  flows) 
is  larger  than  the  other  components.  On  the  other  hand,  when  turbulence  is  created  by  a  grid 
oscillating  in  the  vertical  direction,  the  largest  fluctuations  are  in  the  vertical  direction;  this 
is  not  a  direct  effect  of  the  grid  motion  but  is  due  to  a  selection  process  which  allows  vertical 
motions  to  penetrate  further  into  the  stratified  fluid.  The  mixing  efficiency  correlation  derived 
for  shear  turbulence  does  not  work  in  this  case  and  it  is  necessary  to  create  a  Froude  number 
based  on  the  vertical  component  of  the  turbulence  to  collapse  all  of  the  data.  Turbulent  diffu¬ 
sion  (transport  of  turbulent  energy  through  the  agency  of  the  third  order  correlations)  exhibits 
similar  differences  in  shear-  and  grid-generated  turbulence  and  for  similar  reasons. 

3.2  Oceanic  Upwelling  ^ 

An  interesting  phenomenon  occurs  on  the  eastern  boundaries  of  oceans  i.e.  off  the  west  coasts  of 
continents.  Here,  one  often  finds  stable  high  pressure  areas  (especially  in  local  summer)  which 
drive  long-shore  currents  which,  in  the  northern  hemisphere,  are  from  north  to  south.  The 
Coriolis  force  drives  these  currents  offshore  and  the  surface  water  is  replaced  by  deep  ocean 
water  that  is  colder  and  more  nutrient-laden  than  the  water  it  replaces.  The  deep  water  also  has 
a  different  velocity  than  the  surface  water  and  the  combination  of  velocity  and  density  gradients 
can  lead  to  an  instability  that  makes  the  boundary  between  the  two  fluid  masses  very  irregular. 
This  upwelling  process  is  responsible  for  a  surprisingly  large  part  of  the  mixing  that  occurs  in 
the  ocean;  some  estimates  say  that  half  of  all  mixing  between  surface  and  deep  water  is  due  to 
upwelling. 

Despite  the  fact  that  the  flow  is  turbulent,  many  of  the  properties  of  upwelling  can  be 
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predicted  with  the  aid  of  linear  stability  theory.  It  is  also  possible  to  create  a  laboratory 
experiment  that  mimics  the  process  and  to  perform  simulations  of  it.  Both  of  the  latter  have 
relatively  smooth  surfaces  not  representative  of  the  actual  ocean  but  it  is  possible  to  introduce 
idealizations  of  coastal  geography  and  bathymetry.  These  show  that  capes  and  ridges  cause 
significant  increases  in  the  mixing  that  occurs  in  the  ocean. 
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Numerical  Considerations  in  Large  Eddy 
Simulation  of  Complex  Turbulent  Flows 

Parviz  Moin 
Stanford  University 


We  report  on  numerical  issues  in  large  eddy  simulations  that  have  been 
identified  in  recent  computations  and  theoretical  studies.  These  issues  range 
from  spatial  discretizations  to  the  impact  of  boundary  conditions  and  grid 
resolution.  We  also  describe  the  results  from  several  large  eddy  simulations 
of  complex  turbulent  flows  including  those  performed  with  a  novel  numerical 
technique  based  on  B-splines. 

Computations  of  flow  over  a  circular  cylinder  at  the  sub-critical  Reynolds 
number  of  3900  were  carried  out  with  a  fifth  order  upwind  biased  scheme, 
second  order  central  difference  scheme  and  the  B-spline  method  with  a  zonal 
grid.  One-dimensional  power  spectra  of  the  velocity  fluctuations  in  the  wake 
of  the  cylinder  show  that  the  inherent  dissipation  in  the  upwind  scheme 
leads  to  very  pronounced  damping  of  the  velocity  fluctuations  at  medium 
to  high  frequencies.  The  agreement  of  the  power  spectra  obtained  with  the 
B-spline  method  with  the  data  is  excellent.  The  non-dissipative  second  order 
scheme  agrees  with  the  data  better  than  the  fifth  order  upwind  scheme. 
Near  the  cylinder,  turbulence  statistics  from  the  three  computations  axe  in 
good  agreement  with  each  other,  however,  further  downstream,  the  B-spline 
method  is  in  significantly  better  agreement  with  the  data.  An  interesting 
result  from  these  computations  was  the  demonstration  that  inadequate  grid 
resolution  can  lead  to  early  transition  in  shear  layers  emanating  from  the 
cylinder.  A  cautionary  note  is  that  such  computations  may  show  fortuitous 
agreement  with  experiments  that  also  suffer  from  early  transition  due  to 
external  disturbances. 

Additional  examples  of  large  eddy  simulation  of  complex  flows  that  will  be 
presented  are:  combustion  in  a  co-axial  jet  combustor,  flow  near  the  trailing 
edge  of  a  hydrofoil  and  flow  in  an  asymmetric  diffuser.  In  the  diffuser  flow,  it 
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will  be  demonstrated  that  the  flow  is  very  sensitive  to  inflow  velocity  profile. 
In  the  case  of  the  combustor,  the  exit  boundary  has  a  profound  effect  on 
the  vortex  breakdown  which  is  essential  for  flame  stability  and  mixing  in  the 
combustor. 

Finally,  a  summary  of  some  of  the  current  efforts  aimed  at  making  LES 
an  engineering  tool  will  be  presented.  Specifically,  the  development  of  suit¬ 
able  wall  boundary  conditions  for  LES  and  parallelization  of  LES  codes  will 
be  discussed.  Chapman’s(1979)  estimates  for  the  required  computational  re¬ 
sources  and  the  demonstrated  capabilities  of  DOE’s  ASCI  computers  will  be 
used  to  project  the  near  term  utility  of  LES. 


6 


Insights  for  LES  from  Structure-Based  Turbulence  Modeling 

S.  C.  Kassinos  and  W.  C.  Reynolds 
Department  of  Mechanical  Engineering 
Stanford  University 
Stanford,  CA  94305-3030,  USA 

We  have  spent  the  past  several  years  constructing  improved  one-point  Reynolds- 
Averaged  Navier-Stokes  (RANS)  turbulence  models  based  on  new  structure  concepts 
and  new  one-point  tensors  that  carry  important  structural  information.  From  this 
work  we  have  learned  a  great  deal  about  how  (and  why)  homogeneous  turbulence 
responds  to  different  types  of  mean  deformation.  Insights  from  this  work  are  impor¬ 
tant  in  Large  Eddy  Simulation  (LES)  because,  except  in  the  wall  region,  LES  models 
assume  local  homogeneity. 

The  Smagorinski  model  and  others  like  it  assume  equilibrium  between  the  locally 
homogeneous  sub-grid  turbulent  stress  field  and  the  instantaneous  local  strain  rate  of 
the  resolved  field.  For  some  types  of  deformation,  such  an  equilibrium  might  not  even 
exist.  If  it  does,  the  model  should  take  into  account  the  alteration  of  this  equilibrium 
due  to  the  type  of  mean  deformation,  and  due  to  strong  mean  or  frame  rotation,  but 
these  effects  are  not  properly  included  in  any  model  that  we  have  seen. 

One  key  example  is  the  case  of  turbulence  subjected  to  an  imposed  axisymmetric 
expansion,  say  with  mean  (or  large-eddy)  strain-rates  =  -T  and  S22  —  S33  =  r/2. 
Analysis,  DNS,  and  experiments  all  indicate  that  the  Reynolds  stress  anisotropy  in 
this  flow  is  larger  if  the  deformation  is  slow  than  if  it  is  rapid.  The  anisotropy 
appears  to  grow  continuously  and  shows  no  sign  of  leveling  off  at  an  equilibrium  state. 
Moreover,  upon  removal  of  the  imposed  strain-rate,  the  stress  anisotropy  continues 
to  increase,  rather  than  returning  to  isotropy  as  the  lore  of  turbulence  says  should 
happen.  Since  regions  of  local  momentary  axisymmetric  expansion  occur  all  over 
LES  flows,  it  would  seem  that  sub-grid  models  that  could  capture  this  effect  would 
be  preferable. 

Another  case  where  current  LES  sub-grid  models  are  weak  is  in  flows  with  strong 
frame  or  mean  rotation.  Under  combinations  of  strain  and  rotation  that  lead  to 
elliptical  streamlines  of  the  mean  flow,  the  turbulence  behaves  very  strangely.  The 
turbulence  state  oscillates  and  does  not  equilibrate;  DNS  simulations  (Blaisdell)  shows 
that  the  kinetic  energy  of  the  turbulence  grows,  whereas  quasi-equilibrium  turbulence 
models  of  the  sort  used  in  LES  for  sub-grid  turbulence  predict  decay. 

The  reasons  for  these  seemingly  strange  effects  become  clear  when  one  understands 
our  new  one-point  structure  tensors,  the  way  that  they  axe  related,  and  the  way  that 
they  are  altered  by  mean  deformation.  All  this  should  be  useful  in  considering  the 
way  that  sub-grid  scale  turbulence  should  be  related  to  the  imposed  resolved-scale 
deformation.  The  purpose  of  this  paper  is  to  review  what  has  been  learned  from  our 
structure-based  analyses,  and  to  suggest  how  it  might  be  applied  in  the  development 
of  better  sub-grid  turbulence  models  for  LES. 
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Abstract.  Methods  for  the  direct  and  large-eddy  simulation  of  flows  about  turbo- 
machinery  blading  are  presented.  For  incompressible  flows,  both  direct  and  large- 
eddy  simulation  require  the  solution  of  a  Poisson  equation  for  pressure,  which  is  of¬ 
ten  the  main  source  of  computational  expense  in  simulations  within  complex  bound¬ 
aries.  Marked  computational  advantages  are  described  for  a  flow  which  is  periodic 
in  one  dimension  through  the  use  of  Fourier  techniques.  With  simple  conditions 
fulfilled,  it  is  possible  to  discrete  Fourier  transform  the  Poisson  equation  for  the 
pressure  field,  decomposing  the  problem  into  a  set  of  two-dimensional  problems  of 
similar  type.  We  demonstrate  that  even  when  a  complex  geometry  and  body-fitted 
curvilinear  coordinates  are  used  in  the  other  two  dimensions,  for  instance  to  solve 
flow  around  a  2D  turbine  blade  or  blade  row,  the  resulting  Fourier-transformed  2D 
problems  are  much  more  efficiently  solved  than  the  3D  problem  by  iterative  means. 

We  show  results  for  the  case  of  LES  of  flow  over  an  idealised  turbine  leading 
edge  shape.  The  simulation  shows  laminar  separation  from  the  point  where  the 
surface  curvature  changes,  2D  instabilities  in  the  separated  shear  layer,  break¬ 
down  to  large-scale  3D  structures  prior  to  reattachment.  Of  particular  interest  is 
the  interaction  with  the  wake-like  outer  flow  with  the  re-established  shear  at  the 
wall,  leading  to  a  secondary  boundary-layer  transition  similar  to  bypass  transition. 
The  simulation  is  carried  out  using  general  curvilinear  coordinates  with  the  con- 
tra variant  velocity  components  and  finite  volume  discretisation,  together  with  a 
generalised  dynamic  subgrid  scale  model 


1  Introduction 

Simulations  of  transition  and  turbulence  over  turbomachinery  blades  are  nec¬ 
essarily  three-dimensional.  However  they  frequently  need  to  be  performed  in 
computational  domains  in  which  one  dimension  can  be  assumed  uniform, 
with  a  statistically  homogeneous  flow  solution  in  that  dimension  which  may 
be  assumed  to  be  periodic.  Derivatives  in  the  periodic  direction  can  then  be 
computed  efficiently  through  the  use  of  discrete  Fourier  transforms.  •, 

The  periodic  dimension  is  the  spanwise  dimension  z.  The  requirements  for 
the  use  of  the  Fourier  method  are:  (i)  That  the  unknown  field  whose  solution 
is  sought  may  be  assumed  periodic  in  z,  an  assumption  which  can  be  made  if 
the  turbulent  pressure  field  is  statistically  homogeneous  in  z  and  the  periodic 
length  is  greater  than  twice  the  correlation  length  of  the  field  in  z;  (ii)  that  the 
mesh  used  in  the  z  direction  is  even;  and  (iii)  that  any  differencing  formulae 
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applied  to  the  field  to  generate  the  discretised  equations  are  also  homogeneous 
(this  rules  out  certain  types  of  upwind  differencing  in  z). 

The  gains  in  computational  efficiency  through  the  use  of  Fourier  methods 
are  generally  very  significant.  They  arise  partly  because  of  the  reduction  of 
the  dimension  of  the  problem  and  the  lack  of  any  connection  between  the  so¬ 
lution  for  different  discrete  wave  numbers,  and  also  partly  because  the  higher 
wave-number  problems  have  increased  diagonal  dominance  in  the  solution 
of  Poisson-like  equations.  This  results  in  greatly  accelerated  convergence  of 
the  higher  wave-number  problems  and  a  corresponding  saving  of  computer 
resource. 

A  large-eddy  simulation  of  flow  over  a  NACA  4412  aerofoil  has  been 
performed  by  Kaltenbach  and  Choi  [1]  at  an  angle  of  attack  of  12®  and  a 
chord  Reynolds  number  of  1.64  x  106.  They  used  body-fitted  coordinates 
with  second-order  finite  differences  and  the  contravariant  velocities  discre¬ 
tised  on  a  staggered  mesh  -  a  method  which  is  very  similar  to  that  used 
by  us  for  somewhat  less  ambitious  flows.  Their  mesh  of  638  x  79  x  48  cells 
was  barely  adequate  for  an  LES  of  the  4412  foil,  with  both  streamwise  and 
spanwise  cell  sizes  being  rather  large  in  terms  of  wall  units,  and  a  number 
of  less  satisfactory  results  were  noted.  A  smoothing  filter  was  required  to 
suppress  grid-scale  oscillations  in  the  laminar  flow  upstream,  and  transition 
to  turbulence  was  found  to  occur  immediately  following  the  region  in  which 
the  filter  was  applied.  Values  of  the  shape  factor  H  dropped  more  rapidly 
in  the  simulated  boundary  layer  in  the  mid-chord  region  of  the  aerofoil  than 
is  found  experimentally  [2]  and  near  the  trailing  edge  there  were  substantial 
differences  in  the  turbulence  intensities  between  simulation  and  experiment. 

2  Covariant  Navier-Stokes  equations 

We  are  primarily  interested  in  geometries  that  are  quite  complex  in  two 
dimensions  such  as  a  2D  aerofoil  or  turbine  blade  row.  These  situations  natu¬ 
rally  involve  curved  surfaces  over  which  the  development  of  boundary  layers, 
including  transition,  separation,  reattachment  and  sometimes  relaminarisa- 
tion  is  to  be  followed.  Here  the  natural  approach  is  to  use  general  curvilinear 
coordinates,  and  this  section  describes  one  approach  to  this  aspect  of  the 
simulation.  Note,  however,  that  the  Fourier  techniques  described  below  for 
pressure  solution  also  work  very  well  for  2D  geometries  with  sharp  corners 
such  as  the  flow  over  a  square  cylinder  [3]. 

There  are  two  main  approaches  to  Navier-Stokes  simulation  in  general 
curvilinear  coordinates:  one  using  scalar  generalised  equations  for  the  conser¬ 
vation  of  Cartesian  components  of  mementum;  and  a  second  method,  used 
here,  based  on  the  fully  general  covariant  Navier-Stokes  equations  represent¬ 
ing  the  conservation  of  contravariant  vector  velocity  components. 

The  familiar  NS  equations  become 

&(•/«’)  =  -JgijdjP  -  J( u'V)  j  +  2J(vs'i)  j  (1) 
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where  J  is  the  Jacobian  of  the  coordinate  transformation  and  g,J  is  the  metric 
tensor.  These  geometric  quantities  are  computed  in  advance  of  the  simula¬ 
tion  to  high  accuracy.  The  equations  are  multiplied  by  the  Jacobian  of  the 
transformation  J  in  order  to  cast  the  equations  and  the  discretisation  in 
conservative  form,  p  is  the  pressure  over  density.  The  strain  rate  is  given  by 

sii  =  gikujk+g*ku<k  =  gik(dkuJ  +  rji[ul)  +  i»j,  (2) 

an  expression  that  involves  the  connection  coefficients  which  are  com¬ 
puted  in  advance.  The  symbol  +U+j  indicates  the  addition  of  symmetrising 
terms  identical  to  those  proceeding  but  with  the  indices  i  and  j  interchanged. 
dj  is  the  partied  derivative  with  respect  to  and  the  j  notation  indicates  the 
covariant  derivative,  also  involving  the  connection  coefficients.  The  covariant 
divergences  in  equation  (1),  for  instance,  are 

7(U’V)  J  =  dj(W)  +  ,  (3) 

and  similarly  for  the  viscous  term.  The  mass  conservation  law  in  general 
coordinates  reads 

di  (JV)  =  0,  (4) 

and  hence  one  can  derive  a  generalised  Poisson  equation  for  the  pressure 
simply  by  taking  the  divergence  of  equation  (1). 

The  momentum  advancement  is  explicit  except  for  the  pressure  term 
which  is  solved  by  a  standard  projection  method.  The  spatial  discretisation  is 
performed  by  one  of  the  standard  approaches  used  in  LES,  with  the  geomet¬ 
ric  quantities  computed  to  higher  order  in  space  in  advance  of  the  simulation 
for  several  distinct  points  of  the  staggered  mesh  arrangement.  Note  that  all 
the  geometric  quantities  are  2D  fields,  with  values  independent  of  z,  greatly 
reducing  the  storage  requirements. 

3  Fourier  pressure  solution 

We  simulate  turbomachinery  blades  and  related  flows  using  an  approach  sim¬ 
ilar  broadly  to  that  of  Kaltenbach  and  Choi  [1]  which  allows  a  partial  use  of 
the  direct  Fourier  technique  for  pressure  solution.  Talking  the  divergence  of 
equation  (1),  one  obtains  the  general  version  of  the  familiar  pressure  Poisson 
equation: 

(5) 

rapid 

(6) 


di(Jg'jdjp)  =  R , 

where  R  is  the  divergence  of  the  provisional  mass  flux  field  Ju*. 

The  equation  can  now  be  discrete  Fourier  transformed  in  z  (a  very 
computational  task)  to  obtain  a  set  of  decoupled  equations, 


dx  +  dyJg22  dy 


10 


Peter  R  Voke  and  Zhiyin  Yang 


(This  process  can  be  performed  even  when  the  z.  derivatives  are  replaced  by 
finite-difference  formulae,  provided  z  in  the  simulation  is  periodic  and  has  an 
even  mesh,  though  kz  then  has  a  slightly  different  meaning.) 

The  2D  equations  (6),  one  for  each  value  of  kz,  can  be  solved  very  fast  even 
when  the  geometry  is  complex.  Non-zero  values  of  kz  increase  the  diagonal 
dominance  of  the  system  of  equations,  resulting  in  very  rapid  convergence  of 
any  iterative  scheme  employed  regardless  of  the  type  of  discretisation  used  for 
the  x  and  y  derivatives.  The  small  number  of  low-&*  equations  which  converge 
too  slowly  are  speeded  up  using  acceleration  techniques  such  as  multigrid  [4]. 
By  this  means,  it  is  possible  to  obtain  pressure  solution  times  for  LES  in 
2D  complex  geometries  that  are  almost  as  fast  in  some  cases  (though  not  as 
accurate)  as  the  direct  solvers  available  in  very  simple  geometries. 

Equation  (6)  contains  cross-derivative  terms  which  generally  affect  the 
diagonal  dominance  of  the  matrices  and  slow  the  convergence  of  iterative 
solutions.  These  terms  increase  in  magnitude  with  the  skewness  of  the  mesh 
(they  are  absent  altogether  for. orthogonal  coordinates).  Thus  highly  skewed 
meshes  should  be  avoided  for  the  sake  of  efficiency.  Nevertheless  we  expect 
the  use  of  a  Fourier  technique  to  yield  performance  enhancements  in  the 
pressure  solution  compared  to  a  fully  3D  solver  of  a  least  a  factor  of  5  and 
sometimes  up  to  60. 

4  Dynamic  subgrid-scale  procedure 

A  dynamic  subgrid-scale  procedure  is  employed  in  our  simulations  since  the 
flow  is  Imainra,  transitional  and  finally  fully  turbulent.  A  dynamic  model 
is  a  method  of  accounting  for  the  unresolved  motions  and  their  effect  on  a 
large-eddy  simulation  that  uses  the  information  available  in  the  simulation 
itself  maximally.  We  can  think  of  a  large-eddy  simulation  as  accurately  repre¬ 
senting  filtered  velocity  fields  after  a  large-scale  filter  has  been  applied  to  the 
true  velocity  fields.  When  applied  to  the  Navier-Stokes  equations,  the  filter, 
represented  by  an  overbar,  results  in  equations  that  look  superficially  very 
like  the  time-averaged  (Cartesian)  Navier-Stokes  equations: 

dtlj  __  1  dp  djuiuj)  d  ( duj  duj  \ 

dt  ~  pdxi  dxj  dxj U  \  dxj  dxi) 

The  term  ufuj  cannot  be  computed,  and  the  difference  between  it  and  the 
part  that  can  be  computed,  is  the  subgrid  Reynolds  stress  r<j,  analogous 
to  the  familiar  Reynolds  stress  occuring  in  the  time-average  Navier-Stokes 
equations. 

The  subgrid  Reynolds  stress  represents  the  effects  of  the  subgrid  motions 
on  the  resolved  fields  of  the  LES.  In  contrast  to  the  standard  Reynolds  stress 
whose  length  and  velocity  scales  are  those  of  the  entire  turbulent  flow  field, 
the  length  scale  and  velocity  scale  associated  with  the  subgrid  Reynolds  stress 
can  be  deduced  simply  and  on  a  local  basis  from  the  mesh  and  the  resolved 
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velocity  field  of  the  simulation.  The  length  scale  is  derived  from  the  local 
mesh  {A)  and  the  velocity  scale  is  dictated  by  the  small-scale  motions  on  the 
mesh,  roughly  equivalent  to  the  largest  sub-grid  velocities. 

The  earliest  estimate  for  the  velocity  scale,  based  on  the  local  strain  rate 
scalar  s,  due  to  Smagorinsky  [5],  was  sA.  He  went  on  to  combine  these  scales 
into  the  simple  gradient  diffusion  model: 

T~ij  —  SijTkk/S  —  Vt  Sij  (8) 

v,  =  CA2s  (9) 

s  =  y/ZsijSij  (10) 

C  is  predicted  theoretically  from  the  Kolmogorov  spectrum  in  homogeneous 
isotropic  turbulence  [6]  to  be  0.172. 

This  very  simple  model  proved  surprisingly  successful  for  several  decades. 
The  model  in  fact  has  clear  shortcomings.  Tests  of  its  predictions  against 
directly  simulated  flow  fields  show  that  the  correlation  between  the  tensors 
in  equation  (8)  is  low.  In  shear  flows,  particularly  wall-bounded  flows,  the 
model  is  far  too  strong  and  the  value  of  C  has  to  be  decreased  in  a  pragmatic 
way;  the  reasons  for  this  are  now  quite  well  understood. 

The  dynamic  procedure  is  a  new  approach  in  which  the  value  of  C,  instead 
of  being  adjusted  artificially  based  on  the  experience  of  the  simulator,  is 
deduced  from  the  LES  itself  [7].  This  allows  the  value  of  C  to  respond  to 
local  flow  conditions,  varying  in  time  and  space.  A  base  subgrid-scale  model 
such  as  Smagorinsky  is  needed  first.  We  then  consider  what  would  happen  if 
the  same  LES  flow  field  were  simulated  on  a  coarser  mesh,  effectively  filtering 
the  velocity  field  on  a  sacle  larger  than  that  actually  used  for  the  LES.  This 
coarser  filter,  called  the  test  filter,  is  usually  chosen  to  be  twice  the  size  of 
the  mesh.  The  result  of  applying  the  test  filter  (indicated  by  a  caret)  to  the 
simulated  flow  fields  is  the  twice-filtered  Navier-Stokes  equations: 

d%  1  dp  d  .  .  ( d%  dtj\ 

w=-^-\r+r, I11) 

Once  more  the  filtered  non-linear  term  is  modelled,  but  in  this  case  the  dif¬ 
ference 

Tj  =  Ttfi }  -  1 7i  Uj  (12) 

is  called  the  subtest-scale  stress. 

The  dynamic  procedure  is  based  on  the  assumption  that  the  subgrid  stress 
and  subtest  stress  can  be  modelled  in  a  formally  similar  manner,  by  applying 
the  basis  model  at  both  filter  scales,  A  and  2 A: 

Tij  — "OcAe  / 3  =  CA2  S  Sij  (13) 

Tij  -  SijTkk /3  =  C(2A) 2 1 -tij  .  ( 14) 

Germano  pointed  out  an  exact  identity  relating  the  traceless  part  of  the  test- 
to-grid  transfer  stress,  which  is  computable  in  the  simulation,  to  the  difference 


12 


Peter  R  Yoke  and  Zhiyin  Yang 


between  the  two  modelled  stresses,  the  subtest-scale  and  the  subgrid-scale 
stress: 

(15) 

L\j  =  Lij  -  8ijLkk/ 3  =  Tij  -  SijTkk/ 3  -  (r,-j  -  tu/3)  =  CM{j  ,  (16) 

The  models  for  Tij  and  r,j  from  equations  (13)  and  (14)  are  substituted  into 
(16)  to  define  the  tensor  Mij.  Thus  we  deduce  that  C  must  be 

C  =  L'ij/Mij  ,  (17) 

yielding  five  separate  equations  for  C  since  five  components  of  each  tensor 
on  the  right  hand  side  are  independent  (the  tensors  being  symmetric  and 
traceless). 

The  least-squares  estimate  of  the  optimal  solution  for  C  is  [8] 


^  L'ijMij  > 

<  MkiMki  >  ’ 


(18) 


where  the  angle  brackets  represent  an  average  over  the  homogeneous  direction 
z  in  the  flow  in  which  we  do  not  wish  the  dynamic  procedure  to  generate  local 
variation  of  C.  The  resulting  C  is  a  function  of  time  and  the  inhomogeneous 
coordinates  x  and  y. 

The  dynamic  procedure  has  been  presented  above  in  Cartesian  coordi¬ 
nates  for  the  sake  of  simplicity,  but  its  derivation  can  be  carried  out  in  general 
curvilinear  coordinates  by  exact  analogy  and  without  difficulty.  The  result  is 
that  we  can  deduce  a  dynamic  subgrid  coefficient  C  as 


_  <L'iigikgjiMkl  > 
<  MklgikgjiMkl  >  ’ 


where  L'ij  is  the  traceless  part  of 


xPu?  —  tTuJ  (20) 

and  Mkl  is  similarly  the  contravariant  counterpart  of  Mki  arising  from  the 
difference  between  test-  and  subgrid-scale  modelled  stresses  in  general  coor¬ 
dinates.  A  question  does  arise  regarding  the  variation  of  the  filter  scale  on  a 
curvilinear  mesh,  since  we  treat  the  uniform  mesh  in  the  transformed  plane 
as  defining  the  filter  at  both  grid  and  test  levels.  Clearly  this  does  not  take 
account  in  any  way  of  the  change  of  scale  of  the  mesh  in  physical  space,,  which 
will  affect  the  turbulence  evolution. 

The  practical  computation  of  C  from  equation  (19)  is  considerably  more 
complex  than  in  the  Cartesian  case  (18)  owing  to  the  presence  of  the  met¬ 
ric  tensor  coupling  the  various  components.  A  simple  alternative  that  we 
have  utilised  for  the  simulation  whose  results  are  given  below  is  to  return  to 
equation  (17),  or  rather  its  covariant  equivalent 

C  =  L'ij/Mij ,  (21) 
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and  select  one  of  the  five  equations  to  define  C.  This  is  justified  in  the  case 
of  transition  over  a  leading  edge  where  there  is  a  clear  principal  strain  rate 
and  component  of  the  subgrid  Reynolds  stress,  namely  the  12  component, 

C  —  U 12 /M12  .  (22) 

The  least  squares  form  involving  the  sum  of  all  six  components  of  the  stresses 
is  dominated  by  the  12  component  in  this  case,  so  the  simplification  (22)  leads 
to  almost  the  same  distribution  of  values  of  C.  Of  course  there  is  also  consid¬ 
erable  saving  in  cost  by  using  the  dynamic  procedure  in  this  form,  without 
detectable  change  in  the  quality  of  the  results.  The  dynamic  computation  of 
C  involving  filtering  on  the  test  scale  is  still  quite  expensive  and  is  carried 
out  not  at  every  time  step  but  every  ten  steps,  generating  a  distribution  of 
C  that  is  used  for  the  following  ten  steps. 

In  the  simulation  of  the  transitional  flow  described  below,  the  value  of  C 
increases  from  zero  as  the  flow  becomes  unstable  and  reaches  normal  LES 
levels  when  the  flow  becomes  fully  turbulent.  In  fact  both  Af,J  and  tend 
to  zero  in  laminar  regions,  making  equations  (21)  or  (22)  poorly  conditioned 
or  singular.  The  algorithm  for  computing  C  behaves  unpredictably  in  laminar 
or  nearly  laminar  regions  and  must  be  modified  to  ensure  that  C  approaches 
zero  in  a  sensible  and  controlled  manner.  We  do  this  by  changing  equation 
(22)  to 

C  —  V  12/(Af 12  +  e) .  (23) 

e  is  a  selected  small  quantity  related,  for  instance,  to  the  global  average  of 
M12.  Typically  e  is  of  the  order  of  10-6  of  the  values  of  M 12  occurring  in 
turbulent  regions. 

5  .  Transition  following  a  curved  leading  edge 

We  present  results  from  a  simulation  of  the  flow  over  a  flat  plate  with  a  semi¬ 
circular  leading  edge,  Figure  1.  The  fully  general  coordinate  system  around 
this  geometry  is  needed  since  an  orthogonal  coordinate  system  will  have  an 
inconvenient  junction  between  a  cylindrical  and  a  Cartesian  mesh  region  just 
at  the  critical  point  where  the  curved  leading  edge  joins  the  flat  surface.  This 
is  the  point  at  which  the  laminar  flow  separates.  Any  disturbances  in  this 
region,  including  those  with  a  numerical  origin,  will  affect  the  flow  radically 
and  lead  to  erroneous  conclusions.  A  previous  simulation  by  us  [9]  used  an 
orthogonal  coordinate  system  and  suffered  from  the  problems  indicated,  in¬ 
cluding  rapid  and  excessive  growth  of  instabilities  in  the  separated  shear  layer 
whose  origin  could  be  traced  to  numerical  effects.  This  does  not  occur  in  the 
simulation  with  fully  general  curvilinear  coordinates. 

The  overall  geometry  is  indicated  by  an  instantaneous  flow  field  in  Figure 
1.  The  Reynolds  number  based  on  the  plate  thickness  d  is  3450.  The  simu¬ 
lation  uses  408  (streamwise,  wrapped  round  the  leading  edge)  by  72  (wall- 
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normal)  by  64  (spanwise)  meshes  and  is  performed  using  the  techniques  out¬ 
lined  in  previous  sections,  including  the  generalised  dynamic  subgrid-scale 
model.  The  circular  inflow  boundary  and  the  lateral  boundaries  are  eight 
leading-edge  diameters  (8 d)  distant  from  the  surface,  corresponding  to  a 
blockage  ratio  of  16.  The  inflow  velocity  is  constant,  U0,  and  aligned  with 
the  plate.  The  lateral  boundaries  are  free-slip  but  impermeable.  On  the  out¬ 
flow  boundaries,  9.5tf  downstream  of  the  leading  edge,  we  apply  a  convective 
boundary  condition  based  on  the  mean  streamwise  velocity. 

The  spanwise  dimension  of  the  domain  is  2d.  Some  simulations  have  been 
performed  using  a  spanwise  dimension  of  4 d  without  any  appreciable  change 
in  the  behaviour  of  the  flow.  In  terms  of  wall  units  based  on  the  turbulent 
boundary  layer  downstream  of  reattachment,  the  streawise  mesh  sizes  vary 
from  Ax  =  10  to  30.5,  while  Az  —  9  and  at  the  wall  Ay  =  1  (the  distance 
from  the  wall  to  the  first  specified  value  of  u  and  w.  The  time  step  used  in 
the  simulation  is  0.005t/o/d. 

Statistics  are  gathered  by 'averaging  in  time  once  the  simulation  is  in 
its  fully  developed  state  and  also  over  the  span  direction  and  on  both  sides 
of  the  plate.  The  simulation  is  run  for  4000  time  steps  (200 U0/d)  to  allow 
the  transition  and  turbulent  boundary  layer  to  become  established,  and  the 
statistics  discussed  below  are  then  gathered  over  a  further  30000  time  steps 
(l50Uo/d),  with  a  sample  taken  every  20  steps  (1500  samples). 

Insignificant  small-scale  instability  is  found  around  the  leading  edge,  but 
the  stagnation  point  does  move  from  side  to  side  if  the  division  of  the  mass 
flux  on  either  sides  of  the  plate  is  not  constrained.  The  simulation  can  become 
permanently  asymmetric  with  this  outflow  condition,  so  we  have  constrained 
the  mass  flux  through  the  outflows  to  be  precisely  equal  on  either  side  of  the 
plate  at  every  time  step,  the  total  flux  balancing  the  inflow. 

,  Figures  2  and  3  show  the  mean  and  rms  fluctuating  parts  of  the  streamwise 
velocity  compared  with  experiment  [10]  at  three  streamwise  stations,  one  in 
the  centre  of  the  bubble  ( x/l  =  0.44),  one  just  beyond  the  mean  reattachment 
position  ( x/l  =  1.09  and  one  a  little  further  downstream  (x/l  =  1.64).  Both 
experiment  and  data  have  the  same  low  level  of  free-stream  turbulence  (< 
0.1%)  though  in  the  simulation  this  weak  disturbance  is  imposed  artificially 
through  pseudorandom  forcing  just  upstream  of  the  separation  point.  The 
experimental  blockage  ratio  is  lower  than  that  in  the  simulation,  leading  to 
a  slightly  longer  bubble  length  (2.75 d)  than  in  the  simulation  (about  2.2 d). 
The  experimental  data  are  therefore  compared  at  corresponding  values  of  x/l, 
where  l  is  the  mean  reattachment  length.  For  the  same  reason  all  profiles  are 
plotted  as  functions  of  y/l. 

The  agreement  of  the  mean  profiles  is  excellent,  apart  from  the  velocity  at 
the  edge  of  the  bubble  which  is  higher  in  the  simulation  owing  to  the  different 
blockage  ratio.  However,  in  the  bubble  the  simulation  shows  higher  peaks  of 
u'  occurring  closer  to  the  wall.  The  experimental  data  was  taken  with  a  single 
hot-wire  probe  in  the  shear  layer  and  the  bubble,  and  has  been  discarded  in 
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the  region  where  measured  rms  u'  exceeds  30%  .of  the  local  mean  U.  The 
disagreement  seen  in  this  region  appears  to  arise  from  lack  of  resolution  of 
the  transition  process  in  the  thin  bubble  by  the  LES  rather  than  experimental 
errors,  since  the  position  of  the  peak  fluctuation  relative  to  the  mean  profile  is 
known  to  agree  for  a  wide  variety  of  separation  bubbles  and  is  in  the  position 
found  in  this  experiment  also.  Otherwise  the  agreement  is  excellent. 

Note  that  at  the  x/l  =  1.64  station  the  simulation  shows  a  double  peak 
of  u',  a  new  near-wall  peak  having  appeared  below  the  peak  that  persists 
downstream  of  the  free  shear  layer.  This  effect  is  not  apparent  in  the  ex¬ 
perimental  data  at  this  station.  It  would  indicate  that  two  regions  exist:  a 
wake-like  disturbance  downstream  of  the  bubble  which  is  partially  separated 
from  the  near  wall  production  region;  further  evidence  of  this  double  layer 
structure  and  its  origin  will  be  presented  below. 

Figure  4  shows  the  vertical  ( v ')  fluctuation,  for  which  experimental  data 
are  not  available.  In  figure  5  we  have  the  mean  principal  shear  dU/dy,  to 
illustrate  the  position  of  the  bubble  and  the  rapid  re-establishment  of  high 
near-wall  shear  immediately  after  reattachment,  underneath  the  separated 
shear  layer.  It  is  important  however  to  bear  in  mind  that  these  quantities 
are  time  means  of  highly  fluctuating  instantaneous  flows,  though  the  origin 
of  the  double  layer  structure  is  clear  in  the  two  distinct  mean  shear  regions. 

We  tire  particularly  interested  in  the  product  v'2dU /dy  which  is  the  main 
contribution  to  the  production  of  the  u'v1  stress.  This  is  shown  in  figure  6, 
normalised  by  the  product  of  the  maximum  v 12  (occurring  here  in  the  sep¬ 
arated  shear  layer)  and  the  maximum  dU/dy ,  which,  beyond  reattachment, 
occurs  at  the  wall.  Note  the  slight  evidence  for  the  double  layer  structure  at 
x/l  =  1.64  in  this  figure. 

Figure  7  shows  the  principal  off-diagonal  stress  u'v'.  There  is  no  evidence 
here  for  a  secondary  peak  of  the  stress  close  to  the  wall  even  at  x/l  =  1.64, 
and  in  fact  the  main  contribution  continues  to  be  from  the  separation  bub¬ 
ble  and  its  wake-like  aftermath  to  the  end  of  the  computational  domain. 
However,  looking  at  the  main  contribution  to  the  u  production  u'v' dU/dy 
in  Figure  8,  we  find  a  powerful  near-wall  peak  has  developed  since  the  level 
of  u'v '  stress  from  the  outer  layer  overlaps  in  a  narrow  region  with  the  high 
dU/dy  peaking  at  the  wall.  This  explains  the  origin,  in  statistical  terms,  of 
the  rapid  development  of  a  near  wall  peak  in  u'  beyond  reattachment. 

In  order  to  increase  our  understanding  of  the  processes  involved  in  this 
transition  process,  we  separate  two-dimensional  and  three-dimensional  parts 
of  the  velocity  field  at  a  specific  instant  in  the  simulation.  Figure  9  shows 
the  mesh  in  the  sub-region  we  focus  on,  and  Figure  10  colour  contours  of  the 
streamwise  velocity  component  at  an  arbitrary  z  plane.  The  picture  clearly 
suggests  a  3D  breakdown  of  the  separated  shear  layer  about  two  thirds  of 
the  way  down  the  bubble,  with  the  asymmetry  between  the  two  sides  being 
quite  clear,  suggesting  a  symmetry-breaking  transition  process  is  under  way. 
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Figures  11  to  14  clarify  this.  Figure  11  shows  span-average  streamwise 
velocity  contours  at  the  same  instant,  revealing  that  the  bubble  is  indeed 
two-dimensional  for  the  first  two-thirds  of  its  length  but  that  only  minor 
aspects  of  the  breakdown  process  involve  coherent  2D  processes,  masked  here 
by  the  strong  overall  profile  of  <u>.  The  2D  component  of  the  breakdown 
process  is  more  clearly  apparent  in  Figure  12  which  shows  the  span-average 
vertical  velocity  <v>.  This  picture  indicates  that  weak  rolls  of  fluid  are 
present,  predominant  wallward  motions  being  interspersed  with  regions  of 
motion  away  from  the  wall.  The  structures  do  occur  in  a  region  away  from 
the  wall,  as  if  forming  an  assymmetric  wake  of  the  separtion  region,  and  form 
a  semi-coherent  pattern  with  irregular  streamwise  spacing. 

The  three-dimensional  parts  of  the  velocity  fields  are  given  as  contours 
in  Figure  13  (ti—  <u>)  and  Figure  14  (v—  <v>)).  The  three-dimensional 
disturbances  in  u  are  seen  first  about  half  way  along  the  bubble,  become 
established  by  about  x/l  =  0.7  and  are  highly  significant  at  and  beyond  the 
mean  reattachment  point.  This  confirms  the  general  view  that  3D  breakdown 
promotes  the  reattachment  and  also  suggests  that  the  highly  variable  position 
of  reattchment  in  time  and  z  is  related  to  the  chapotic  nature  of  the  3D 
breakdown  locally.  Figure  14  reveals  that  vertical  3D  motions  occur  later, 
being  first  seen  just  upstream  of  mean  reattachment,  and  are  less  pronounced, 
with  an  irregular  streamwise  structure.  This  is  also  true  for  span-wise  motions 
(not  shown).  The  3D  motions  however  approach  the  wall  more  closely  than 
the  2D  rolls,  and  the  u  and  v  pictures  are  related,  there  being  a  strong 
suggestion  of  the  correlation  between  negative  u  (yellow-red  in  Figure  13)  and 
motion  away  from  the  wall  (blue  above  the  plate  and  red  below  in  Figue  14). 
The  ingress  of  of  these  3D  vertical  fluctuations  into  the  very  near  wall  region 
where  the  instantaneous  shear  du/dy  is  also  high  supports  our  view  of  the 
mechanism  of  transition. 

We  have  presented  strong  evidence  for  a  simple  conceptual  model  of  the 
process  of  transition  in  this  type  of  flow  from  these  data.  The  separation 
is  two-dimensional  and  laminar.  Three-dimensional  breakdown  occurs  in  the 
second  half  of  the  bubble,  with  the  growth  of  3D  vertical  motions  resulting  in 
turbulent  transport  leading  to  reattachment  which  is  highly  variable  in  time 
and  streamwise  and  spanwise  position.  The  turbulent  second  half  of  the  bub¬ 
ble  gives  rise  to  a  qusi-regular  asymmetric  wake-like  disturbance  propogating 
downstream,  clearly  separate  from  the  wall  layer.  Three-dimensaional  distur¬ 
bances  within  this  ‘wake’,  however,  approach  closer  to  the  wall  and  impinge 
on  the  growing  shear  of  the  near  wall  layer,  promoting  a  second  transition 
process  that  is  similar  in  most  respects  to  bypass  transition  stimulated  by 
similar  levels  of  free-stream  turbulence  [11]. 

It  is  the  relatively  small  scale  of  these  3D  incursions  that  allows  them  to 
reach  in  close  to  the  wall  and  interect  with  the  high  shear  layer,  lifting  low 
speed  fluid  away  from  the  wall  and  pushing  higher  speed  fluid  towards  it.  The 
creation  of  irregularly  spaced  streaks  (not  shown  here)  in  the  near  wall  region 
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is  thought  to  lead  to  redistribution  of  fluctuations  through  pressure-strain 
interaction,  producing  vertical  (v)  fluctuations  and  leading  to  self-sustained 
turbulence  in  the  boundary  layer. 

6  Discussion 

Wew  have  presented  a  numerical  simulation  of  the  transition  to  turbulence 
occurring  through  separation  bubbles  on  a  plate  with  a  semi-circular  leading 
edge.  The  boundary  layers  develop  in  the  densely  meshed  regions  on  either 
side  of  the  plate,  and  are  fully  turbulent  well  before  they  reach  the  two 
outflow  boundaries.  The  layers  on  either  side  of  the  plate  are  statistically 
independent,  just  as  they  would  be  in  the  more  complex  case  of  a  turbine 
blade.  The  flow  separates  and  transition  takes  place  as  the  separated  flow 
becomes  unstable,  first  in  a  two-dimensional  mode  but  passing  into  a  three- 
dimensional  breakdown  which  results  in  rapid  reattachment. 

The  large-scale  instabilities'  originating  in  the  bubble  subsequently  dis¬ 
turb  the  reattached  boundary  layer.  As  the  mean  shear  at  the  wall  is  re¬ 
established,  large-scale  disturbances  in  the  outer  region  of  the  boundary  layer 
coming  downstream  rather  like  a  wake  from  the  separation  bubble  stimulate 
a  secondary  transition  near  the  wall.  The  mechanisms  at  work  here  appear 
to  be  similar  to  those  that  operate  during  bypass  transition  in  the  presence 
of  free-stream  turbulence  [11],  [9]. 

7  Conclusions 

Methods  for  performing  large-eddy  and  direct  numerical  simulations  using 
general  curvilinear  coordinates  have  been  described,  together  with  the  use 
of  Fourier  techniques  for  efficient  pressure  solution  in  this  context.  We  have 
shown  how  it  is  possible  to  combine  iterative  methods  with  discrete  Fourier 
transforms  to  obtain  very  acceptable  pressure  solution  times,  comparable 
with  those  obtained  by  direct  methods  in  simpler  geometries.  Results  for  a 
transitional  flow  using  general  coordinates  in  two  dimensions  demonstrate 
the  efficacy  of  the  approach.  These  methods  are"  currently  in  use  for  the 
simulation  of  turbine  blade  flows  and  will  also  be  used  in  the  LES  of  aerofoils 
in  future. 
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Figure  Captions 

Figurel.  Instantaneous  flow  field  about  the  plate  with  curved  leading  edge: 
whole  domain.  Note  the  fine  mesh  in  the  region  of  primary  interest  close  to 
the  surface. 

Figure  2.  Mean  flow  velocity  U  at  three  positions  measured  from  the  blend 
point.  Left  to  right,  x/l  =  0.44, 1.09, 1.64. 

Figure  3.  RMS  streamwise  velocity  fluctuation  u'  at  x/l  =  0.44, 1.09, 1.64. 

Figure  4.  RMS  wall-normal  velocity  fluctuation  vf  at  *//  =  0.44, 1.09, 1.64. 

Figure  5.  Mean  principal  shear  dU/dy  at  x/l  =  0.44, 1.09, 1.64. 

Figure  6.  Stress  production  term  v'2  dU/dy,  normalised  by  peak  v'2  and  peak 
dU/dy ,  at  x/l  =  0.44, 1.09, 1.64. 

Figure  7.  Principal  Reynolds  stress  uV  at  x/l  =  0.44, 1.09, 1.64. 

Figure  8.  Fluctuation  production  term  u'v1  dU/dy ,  normalised  by  peak  u'v' 
and  peak  dU/dy,  at  x/l  =  0.44, 1.09, 1.64. 

Figure  9.  Mesh  in  the  subdomain  depicted  in  figures  10  to  14. 

Figure  10.  Contours  of  streamwise  velocity  u  at  a  single  instant  for  a  se¬ 
lected  spanwise  plane.  Blue,  positive;  red,  negative. 

Figure  11.  Contours  of  span-averaged  streamwise  velocity  <u>  at  a  single 
instant.  Blue,  positive;  red,  negative. 

Figure  12.  Contours  of  instantaneous  span-averaged  vertical  velocity  <t>>. 
Blue,  positive;  red,  negative. 

Figure  13.  Contours  of  3-dimensional  part  of  the  streamwise  velocity  field 
u—  <u>  at  one  instant  for  a  selected  spanwise  plane.  Blue,  positive;  red, 
negative. 

Figure  14.  Contours  of  the  3-dimensional  part  of  the  vertical  velocity  field 
v~  <v>  at  one  instant  for  a  selected  spanwise  plane.  Blue,  positive;  red, 
negative. 
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Istanbul  1998:  Voke  &  Yang,  Figure  1 


Istanbul  1 998:  Voke  &  Yang,  Figure  2 
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Istanbul  1998:  Voke  &  Yang,  Figure  3 


Istanbul  1998:  Voke  &  Yang,  Figure  4 


Istanbul  1998:  Voke  &  Yang,  Figure  5 
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Istanbul  1 998:  Voke  &  Yang,  Figure  9 


24 


Istanbul  1998:  Voke  &  Yang,  Figure  10 
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Istanbul  1998:  Voke  &  Yang,  Figure  11 
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Istanbul  1998:  Voke  &  Yang,  Figure  12 
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Istanbul  1998:  Voke  &  Yang,  Figure  13 


1.5 

1 


0.5 

0 

-0.5 


-1 

-1.5 


-1 


-0.5 


0 


0.5 


1 


Istanbul  1998:  Voke  &  Yang,  Figure  14 
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Direct  Numerical  Simulation  (DNS)  is  the  most  ac¬ 
curate,  but  also  the  most  expensive,  way  of  comput¬ 
ing  complex  turbulent  flows.  In  view  of  the  com¬ 
putational  complexity  of  DNS,  our  first  concern  is 
to  reduce  the  computational  costs  as  far  as  we  can 
get.  This  implies  -  among  others  -  that  the  number 
of  grid  points  has  to  be  kept  as  small  as  possible. 
To  use  the  lowest  possible  number  of  grid  points, 
spatial  discretization  methods  for  the  Navier-Stokes 
equations  need  to  be  strained  to  their  limit. 

On  nonuniform  grids  various  ways  exist  to  discretize 
convective  and  diffusive  operators.  In  [1],  we  have 
proposed  to  apply  a  high-order  finite-volume  dis¬ 
cretization  method  that  preserves  the  spectral  prop¬ 
erties  of  the  convective  and  diffusive  operators,  i.e. 
convection  <->  skew-symmetric;  diffusion  *■>  symmet¬ 
ric  positive  definite.  It  is  our  experience  that  in  this 
way  the  error  of  the  convective  discretization  does 
not  interfere  with  the  diffusion  on  the  smallest  length 
scales.  The  energy  of  the  discrete  system  is  con¬ 
served,  if  the  physical  dissipation  is  turned  off.  The 
time-advancement  is  carried  out  by  an  explicit  one- 
leg  method.  For  more  details  about  the  numerics,  we 
refer  to  [1]. 

We  have  applied  this  numerical  method  to  simulate 
a  number  of  complex  flows.  The  problem  considered 
in  this  paper  is  the  fully  developed  flow  in  a  channel, 
where  a  matrix  of  cubes  is  placed  regularly  at  one 
wall  of  the  channel.  The  matrix  consists  of  25  x  10 
cubes  in  the  streamwise  and  the  spanwise  direction, 
respectively.  The  width  of  the  channel  is  ZAh,  where 
h  denotes  the  width  of  a  cube.  The  pitch  of  the 
cubes  is  4 h,  in  both  the  streamwise  and  the  spanwise 
direction.  The  Reynolds  number  (based  on  the  width 
of  the  channel  and  the  bulk  velocity)  is  equal  to  Re  = 
13,000. 

Flow  measurements  by  Meinders  &  Hanjalic  [2]  have 
showed  that  the  influence  of  the  in-  and  outlet  can 
be  neglected  around  the  18th  row  of  cubes  (counted 
from  the  inlet).  Hence,  to  simulate  the  flow  there, 
we  may  confine  the  computational  domain  to  a  sub¬ 
channel  unit  of  dimension  4.hxZAhx4h  with  periodic 
boundary  conditions  in  the  stream-  and  spanwise  di¬ 
rection.  Figure  1  displays  a  sub-channel  unit. 

The  flow  through  the  sub-channel  was  one  of  the 
test  cases  at  the  6th  ERCOFTAC /I AHR/COST 
Workshop  on  Refined  Flow  Modeling  [3].  Four 
groups  have  presented  the  results  of  their  Reynolds- 
average  Navier-Stokes  (RANS)  computations.  Both 


Figure  1:  Top- and  side-view  of  a  sub-channel  unit.  Shown 
is  an  instantaneous  flow  field  (of  the  1003  simulation)  at  two 
planes  through  the  centre  of  the  cubes. 


high-  and  low-Re-number  RANS-models  were  ap¬ 
plied.  The  RANS  computation  that  agreed  the  best 
with  the  available  experimental  data  used  a  67  by  72 
by  57  grid.  We  have  computed  solutions  of  the  in¬ 
compressible  Navier-Stokes  equations  (without  using 
any  turbulence  model)  on  a  number  of  grids.  Fig¬ 
ure  2  shows  a  comparison  of  the  mean  streamwise 
velocity  of  our  603  and  the  1003  Navier-Stokes  com¬ 
putations,  that  of  the  best  RANS  computation  and 
the  experimental  data  of  Meinders  &  Hanjalic.  On 
a  corresponding  grid,  the  mean  velocities  computed 
from  the  Reynolds  averaged  Navier-Stokes  equations 
agree  less  with  the  experimental  data  than  the  re¬ 
sults  of  the  603  Navier-Stokes  computation.  The  ve¬ 
locity  profiles  of  the  RANS  computation  are  much 
too  smooth.  In  addition,  the  maxima  of  the  veloci- 
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Figure  2:  Comparison  of  the  mean  streamwise  velocity  at 
half  cube  height.  The  horizontal  corresponds  to  the  span- 
wise  direction.  The  dashed  vertical  lines  are  lines  of  sym¬ 
metry.  The  geometry  is  drawn  to  scale.  The  velocity  scale 
is  shown  for  the  uppermost  profiles  only;  the  other  profiles 
have  the  same  vertical  scale. 

ties  are  located  in  the  symmetry-plane  between  two 
cubes,  which  is  in  distinct  disagreement  with  the  ex¬ 
perimented  data. 

First-  and  second-order  statistics  obtained  from  the 
1003  simulation  at  the  cross  section  of  the  channel 
that  bisects  a  cube  are  compared  to  the  available 
experimental  data  in  Figure  3.  Here,  the  averages 
are  computed  over  200  (non-dimensional)  time  units. 
The  profiles  of  the  mean  streamwise  velocity  and  the 
mean-square  of  the  fluctuating  streamwise  velocity 
are  in  good  agreement  with  the  experiments,  except 
in  the  front  of  a  cube,  where  some  discrepancies  be¬ 
tween  the  mean-squares  of  the  fluctuating  stream- 
wise  velocities  exist. 

So,  in  conclusion,  the  1003  simulation  reproduces 
the  turbulent  fluctuations  reasonably  well,  whereas 
the  603  simulation  displays  some  shortcomings.  This 
naturally  poses  the  question:  how  maiiy  grid  points 
are  required  for  a  numerical  simulation  of  the  flow 
under  consideration,  without  any  turbulence  model? 
To  answer  this  question,  we  have  also  performed  a 
1443  simulation.  With  the  help  of  the  solutions  on 
three  different  grids,  GO3,  1003  and  1443,  the  con¬ 
vergence  of  the  numerical  approximation  upon  grid 
refinement  is  addressed.  Here,  the  experimental  re¬ 
sults  of  Meinders  k  Hanjalic  [2]  form  the  frame  of 
reference.  The  spectra  of  the  603,  1003  and  1443 
simulations  are  compared,  to  identify  the  scales  of 
motion  that  are  significant  for  the  first-  and  second- 
order  statistics.  In  addition,  near  the  flat  wall  of  the 
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Figure  3:  A  comparison  of  first-order  statistics  (upper  pic¬ 
ture)  and  second-order  statistics  (lower  picture)  of  the  DNS 
with  experimental  data.  Shown  are  the  mean  streamwise 
velocity  u  (upper  picture)  and  u'u'  (lower  picture)  in  the 
plane  parallel  to  the  streamwise  direction  that  bisects  the 
cubes.  The  continuous  lines  correspond  to  the  DNS;  the 
experimental  data  is  depicted  by  the  dots. 

channel,  the  numerical  solutions  are  also  compared 
to  those  of  the  DNS  of  Kim  et  al.  [4]  and  those  of 
Gilbert  &  Kleiser  [5].  Up  to  y+  ~  250,  the  1003  and 
1443  solutions  obey  the  log  law  accurately,  where  we 
have  taken  5.5  for  the  additive  constant  and  2.5  for 
the  coefficient. 

Finally,  the  budgets  for  the  Reynolds  stresses  are 
computed  from  flow  fields  of  the  1443  simulation. 

References 

1.  R.W.C.P.  Verstappen  and  A.E.P.  Veldman,  J.  En- 
gng.  Math.  32  (1997)  143-159. 

2.  E.M.  Meinders,  K.  Hanjalic,  The  flow  in  a  matrix  of 
cubical  protrusions  placed  in  a  fully  developed  low 
Re  number  channel  flow,  submitted. 

3.  K.  Hanjalic  and  S.  Obi  (eds.)  Proc.  6th  ER- 
COFTAC/IAHR/COST  Workshop  on  Refined  Flow 
Modeling,  Delft  University  of  Technology  (1997). 

4.  J.  Kim,  P.  Moin  and  R.  Moser,  J.  Fluid  Mtch.  177 
(1987)  133-166. 

5.  N.  Gilbert  and  L.  Kleiser,  Proc.  Turbulent  Shear 
Flows  8  (1991). 


31 


International  Workshop  on 

INDUSTRIAL  AND  ENVIRONMENTAL  APPLICATIONS  OF  DNS  AND  LES,  Istanbul  1998 


Large  eddy  simulation  of  a  transitional  flow 
over  a  backward  facing  step  with  boundary  layer  manipulations 

G.  Barwolff,  TU  Berlin,  FB  Mathematik 
Germany 


We  consider  a  backward  facing  step  channel  and  we  investigate  the  flow  for  a  Reynolds  number  of 
3000  built  with  the  velocity  u0  of  the  block  inflow  profil  and  the  step  height  H  (see  fig.  1).  The  top 
of  the  channel  is  considered  as  a  slip  wall  and  in  the  lateral  direction  a  periodic  behavior  is  assumed. 
With  the  aim  of  the  drag  reduction  or  a  decreasing  reattachment  lenght  in  the  wake  of  the  step 
acoustic  manipulations  of  the  boundary  layer  in  front  of  the  step  were  performed  by  experiments  and 
numerical  simulations  ([1],  [2]  and  [5]).  It  was  found  a  good  agreement  of  the  experimental  and  nu¬ 
merical  results,  especially  in  the  case  of  the  reattachment  length  but  also  related  to  the  mean  velocity 
field  and  the  rms-profiles.  The  numerical  investigations  are  done  as  direct  numerical  simulations  and 
large  eddy  simulations.  For  the  direct  numerical  simulations  about  10  million  grid  points  are  needed 
for  the  spatial  discretization.  In  the  case  of  the  large  eddy  simulations  it  is  possible  to  work  with 
about  500  thousands  of  spatial  grid  points.  We  use  a  subgrid  scale  model  of  Germano  type  following 
Akselvoll/Moin  [3].  During  the  numerical  simulations  the  sufficiency  of  LES  to  describe  the  flow  over 
the  backward  facing  step  physical  correctly  could  be  shown. 

Beside  the  accoustic  boundary  layer  manipulation  over  slits  with  certain  frequencies  experiments  are 
done  with  oblique  backward  facing  steps  or  with  moving  boundaries  to  simulate  oblique  geometries. 
Based  on  the  experiences  of  the  above  discussed  numerical  simulation  of  a  straight  backward  facing 
flow  large  eddy  simulations  with  moving  boundaries  are  under  consideration. 

We  investigate  both  a  moving  upstream  boundary  in  front  of  the  step  (lateral  velocity  vu)  and  a 
moving  boundary  behind  the  step  in  the  downstream  region  of  the  bottom  of  the  channel  (lateral 
velocity  Vd,  see  fig.  2).  The  simulations  are  done  for  a  wide  range  of  lateral  boundary  velocities,  i.e. 
from  viatcrai  —  0.5  u0  to  viaterai  —  2.0  iiq  with  uq  a s  the  inflow  profile  velocity. 


Figure  1:  channel  situation  Figure  2:  Tu  and  Vd  as  a  moving  boundaries 


In  fig.  3  the  mean  wall  shear  stress  velocity  uT  of  the  case  vu  =  2  u0  and  vd  =  0  is  shown  compared 
to  those  of  the  neutral  case  vu  =  0  uq  and  Vd  =  0  in  fig.  4. 


Figure  3:  uT  of  the  emulated  oblique  step  flow 


32 


The  numerical  simulations  are  done  with  a  finite  volume  method  on  staggered  grids  which  is  paral¬ 
lelized  on  a  Cray  T3D/T3E  using  fast  Cray-specific  shared  memory  transfer  utilities  or  the  platform 
independend  MPI  library  ([4]).  The  method  is  of  second  order  in  space  and  time.  The  mass  conser¬ 
vation  per  time  step  was  realized  by  a  pressure-velocity  iteration  method. 

The  parallelized  numerical  model  (DNS/LES)  was  validated  by  comparisons  of  the  numerical  DNS/LES 
results  of  [5]  and  the  experimental  data  of  [1]  and  [6].  The  following  figures  5-8  show  the  result  of  com¬ 
parison  for  the  mean  velocities  and  the  rms-values.  For  these  comparisons  a  neutral,  non-manipulated 
backward  facing  step  flow  for  the  transitional  Reynolds  number  3000  was  considered.  The  experimen¬ 
tal  results  were  produced  by  LDA-measurements. 


Figure  5:  umean  at  i  =  5 H 


Figure  6:  vmean  at  x  =  5 H 


Figure  7:  urm,  at  x  =  5 H 


Figure  8:  urm,  at  x  =  5 H 


The  performance  of  the  parallelized  code  (1.3  *  107  gridpoints,  100  timesteps,  25  pressure-velocity 
iterations)  on  different  mpp  systems  is  shown  in  the  table  1. 


system 

#procs 

mflops 

tC90/< 

J90 

1 

16428 

92 

0.27 

1.00 

J90 

4 

4173 

357 

1.03 

3.94 

J90 

8 

2244 

673 

1.95 

7.32 

J90 

16 

1390 

1087 

3.15 

11.81 

C90 

4380 

338 

1.00 

system 

#procs 

proc*  *  procy 

s?  iffi 

mflops 

tC9o/4 

Sp  =  tT3D(32)/t 

32 

8*4 

3329 

443 

1.31 

1.00 

T3D 

64 

8*8 

1671 

886 

2.62 

1.99 

T3D 

128 

16  *  8 

871 

5.03 

3.82 

T3D 

256 

16  *16 

437 

3387 

10.02 

7.61 

Table  1:  Performance  of  the  FV  code  on  different  mpp  systems 


The  results  of  the  reattachment  lenght  reduction  in  the  case  of  the  moving  upstream  boundary  in 
front  of  the  step  show  a, good  agreement  with  the  known  experimental  measurements.  The  large  eddy 
simulations  with  the  moving  boundary  behind  the  step  are  in  the  beginning  and  the  first  analysis 
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shows  no  remarkable  influences  of  the  lateral  boundary  velocity  v&  to  the  position  of  the  reattachment 
point  Xr  •  A  significant  reduction  of  the  reattachment  Ienght  was  not  found.  The  comparison  of  the 
numerical  results  and  the  experiments  of  such  a  channel  flow  with  a  moving  channel  bottom  behind 
the  step,  which  are  done  now,  will  be  discussed. 
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People  became  concerned  about  a  possible  impact  of  aircraft  emissions  upon  the  atmosphere 
and  the  climate.  The  concern  is  based  on  the  facts  that  the  world-wide  commercial  air  traffic 
grows  strongly  (and  is  predicted  to  grow  strongly  in  the  next  decade)  and  it  is  the  only  direct 
anthropogenic  emission  source  at  cruising  altitudes  where  the  background  concentration  of  the 
trace  gases  is  low  and,  hence,  an  efficacy  is  possibly  large.  To  assess  the  impact  of  aircraft 
emissions  upon  the  atmospheric  state  at  a  global  scale  it  is  prerequisite  to  know  how  the 
emissions  get  mixed  and  dispersed  in  the  immediate  wake  of  the  aircraft,  since  chemical  species 
transformations  and  physical  phase  changes  of  the  exhaust  often  evolve  nonlinearly  and,  thus, 
may  depend  heavily  on  the  mixing  properties  in  the  wake.  Mixing,  in  turn,  is  controlled  by 
the  dynamics  of  the  wake  flow.  The  objective  of  the  work,  therefore,  is  to  provide  a  detailed, 
quantitative,  and  validated  description  of  the  wake  dynamics  and  its  effect  on  mixing  and 
diffusion  of  the  exhaust  from  the  scale  of  the  airplane  to  the  range  of  atmospheric  mesoscale 
flow.  Large-eddy  simulations  are  performed  to  calculate  the  wake  flow  embedded  m  a  stably 
stratified  and  turbulent  atmosphere.  The  numerical  simulation  data  will  be  compared  to  data 
collected  in-situ  and  remotely  during  various  measurement  campaigns. 

The  prime  result  is  that  the  atmospheric  eddies  disturb  the  parallel  vortex  tubes  by  advection. 
Once  distorted  at  that  scale,  the  vortex  induced  velocities  amplify  the  disturbance  according  to 
Crow  s  instability.  The  vortex  tubes  link  after  about  1.5  minutes  and  fond  rings;  the  continuous 
trail  of  exhaust  is  reorganized  in  a  row  of  single  puffs.  Without  any  atmospheric  turbulence,  the 
wake  does  not  experience  a  sinusoidal  instability  but  the  parallel  vortex  tubes  approach  and 
start  to  dissolve  after  2  minutes  when  they  touch.  The  dissolution  is  triggered  by  the  small- 
scale  turbulent  friction  owing  to  the  boundary-layer  turbulence  of  the  aircraft.  The  exhaust  trail 
remains  aligned  along  the  flight  track.  For  both  calm  and  turbulent  situations  the  entrainment 
rate  drops  in  the  interval  between  20  s  and  100  s  by  2  to  2.5  orders  of  magnitude  and  increases 
again  when  the  vortices  collapse. 
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Another  problem,  rises  in  view  of  the  expected  growth  of  air  traffic:  Increasing  demands  on 
the  capacity  and  safety  of  large  airports  have  to  be  faced,  since  aircraft  wake  vortices  (WV) 
may  exert  a  serious  danger  on  following  aircraft  if  the  separation  between  leading  and  following 
aircraft  is  not  sufficient.  Among  others  these  demands  require  new  ways  to  develop  dynamic 
spacing  criteria  for  approach  and  landing.  Since  at  many  airports  aircraft  have  to  join  the  glide 
slope  already  several  miles  before  landing,  the  area  which  has  to  be  controlled  is  quite  large. 
That  is,  we  need  to  know  where  and  how  long  WV  can  exist  along  the  glide  path.  Preferably, 
we  would  like  to  identify  the  atmospheric  conditions  for  which  WV  do  not  exert  any  danger 
to  following  aircraft.  Therefore  the  questions  whether  and  under  which  atmospheric  conditions 
WV  can  stop  their  descent  or  even  rise  again  are  of  primary  interest. 

A  wake  vortex  warning  system  (WVWS)  was  developed  for  the  airport  in  Frankfurt  in  order  to 
run  the  closely  spaced  parallel  runways  separately  at  appropriate  meteorological  conditions.  The 
WVWS  predicts  the  propagation  and  lifespan  of  wake  vortices  in  a  safety  area  of  80m  height 
based  on  statistical  methods.  This  particular  height  of  80m  was  chosen  since  measurement 
campaigns  had  shown  that  WS  do  not  reach  this  height  due  to  the  rebound  caused  by  the 
interaction  of  wake  vortices  with  the  ground.  Pilot  associations  recently  argued  that  also  the 
updrafts  of  thermals  in  a  convective  atmospheric  boundary  layer  (CBL)  could  cause  WS  to  stall 
or  even  to  rise  again  to  the  glide  slope  during  the  final  approach.  We  have  performed  LES  of 
the  evolving  CBL  where  the  wake  vortices  of  an  aircraft  were  superimposed  on  the  windfield. 
The  results  indicate  that  WV  actually  may  rise  in  the  CBL  when  the  velocity  of  the  updraft 
exceeds  the  induced  descent  speed  of  the  WV.  But  at  the  same  time  they  axe  strongly  deformed 
by  the  large-scale  velocity  field  of  the  CBL.  Moreover,  the  decay  of  the  WV  is  considerably 
accelerated  due  to  the  increased  turbulence  levels  in  the  thermals.  To  quantify  these  qualitative 
results  a  comparison  of  the  rolling  moment  being  exerted  on  an  aircraft  crossing  the  domain 
on  various  paths  will  be  presented  for  the  CBL,  a  weakly  turbulent,  and  a  stable  atmosphere. 


Usually,  the  pair  of  aircraft  wake  vortices  sinks  by  mutual  velocity  induction.  However,  close  to 
the  ground  they  may  rise  again  due  to  the  interaction  of  one  vortex  with  the  friction  layer  of 
the  bottom.  It  has  recently  be  found  that  layers  of  wind  shear  may  have  a.. similar  effect.  Again, 
these  rising  vortices  may  be  hazardous  to  following  aircraft.  We  will  furthermore  present  LES 
results  on  the  impact  of  shear  layers  (e.g.  low-level  jets)  upon  the  bouncing  properties  of  wake 
vortices  when  they  pass  through  such  a  layer. 
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M.  Paraschiovoiu  (U.  Toronto),  R.  M.  Lewis  (ICASE) 

In  the  simulation-based  design  process  the  engineer  must  make  numerous 
appeals  to  the  underlying  numerical  simulation  in  order  to  calculate  the 
outputs  —  the  system  performance  metrics  of  interest,  such  as  li  ,  rag, 
or  contaminant  dispersion  -  for  different  values  of  the  design  vanabte. 
The  requirements  on  the  numerical  approximation  are  thus  twofold:  the 
approximation  must  be  sufficiently  coarse  —  and  hence  inexpensive  so  as 
to  permit  repeated  evaluation;  the  approximation  must  be  sufficiently  fine 
so  that  the  numerical  prediction  of  the  desired  outputs  is  representative  of 
the  true  performance  of  the  system. 

A  posteriori  error  control  offers  great  promise  in  reconciling  these  often 
conflicting  requirements.  A  posteriori  analysis  is  composed  of  two  critical 
ingredients:  an  estimation  procedure  which  inexpensively  assesses  the  error 
in  a  particular  numerical  solution;  and  an  adaptive  refinement  proce  ure 
which  exploits  this  error  information  to  optimally  improve  the  numerical 
solution.  The  objectives  of  a  posteriori  error  control  are  twofold:  to  elim¬ 
inate  numerical  uncertainty  —  arguably  the  single  largest  impediment  to 
widespread  adoption  of  simulation-based  design;  and  to  improve  numerical 
efficiency  —  thus  permitting  much  more  extensive  design  exp  oration,  n 
fact,  greater  certainty  is  a  prerequisite  for  greater  efficiency,  we  may  con 
sider  a  less  expensive  (or  even  the  least  expensive)  discretization  only  it  the 
associated  error  can  be  quantified,  and  hence  constrained  and  controlled;  we 
are  no  longer  compelled  to  choose  either  certainty  or  efficiency  ot  can 
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be  achieved. 

Unfortunately,  in  all  earlier  a  posteriori  error  analysis  techniques,  either 

—  in  implicit  approaches  [5,2,1]  —  the  measure  of  the  error  is  not  related 
to  the  actual  engineering  outputs  of  interest,  or  —  in  explicit  approaches 
[3,4,16]  —  the  error  estimates  for  the  engineering  outputs  of  interest  in¬ 
volve  numerous  undetermined  or  uncertain  constants  or  functions ;  in  both 
cases,  quantitative  confirmation  —  and  hence  both  certainty  and  efficiency 

—  is  seriously  compromised,  and  the  relevance  to  engineering  design  greatly 
reduced. 

The  a  posteriori  error  control  method  that  we  are  developing  [10,12,13,8] 
is  certainly  indebted  to  these  earlier  techniques  —  in  particular,  implicit 
Neumann-subproblem  approaches  [5,2,1]  —  for  several  important  concep¬ 
tual  and  mathematical  ingredients,  especially  duality  and  hybridization. 
However,  our  method  considerably  generalizes  these  techniques,  thereby 
providing  a  new  and  critical  “enabling  technology”:  the  ability  to  obtain 
inexpensive,  sharp,  rigorous,  and  quantitative  (“constant-free”)  bounds  for 
the  numerical  error  in  the  engineering  outputs  of  interest.  Our  method  is 
thus  directly  relevant  to  the  design  process,  and  should  lay  the  foundation 
for  systemic  application  of  a  posteriori  error  control  within  the  engineering 
context. 

Our  initial  formulation  focused  on  symmetric  coercive  problems  (e.g., 
Poisson  and  Linear  Elasticity  [10,12,13]),  and  nonsymmetric  coercive  prob¬ 
lems  (e.g.,  Convection-Diffusion  [10,12,8]),  as  well  as  certain  constrained 
problems  (the  Stokes  equations  [11],  central  to  hydrodynamics).  However, 
we  have  recently  proposed  a  more  general  formulation  that  greatly  expands 
the  types  of  equations  and  outputs  amenable  to  our  approach.  In  partic¬ 
ular,  both  noncoercive  and  nonlinear  elliptic  [9]  (and  stabilized  hyperbolic 
[7])  problems  can  now  be  addressed,  with  only  a  minor  loss  in  certainty  rela¬ 
tive  to  the  coercive  case.  More  precisely,  for  noncoercive  equations  we  must 
require  an  additional,  rather  weak,  hypothesis  related  to  the  relative  con¬ 
vergence  rates  and  magnitudes  of  the  L2  and  Hl  errors;  strictly  speaking, 
we  therefore  provide  only  asymptotic  bounds,  although  in  practice  uniform 
bounds  are  almost  always  obtained. 

We  have  now  demonstrated,  both  theoretically  and  empirically,  that  our 
bound  (and  associated  adaptive  refinement)  procedure  —  rigorous,  quan¬ 
titative,  and  directly  relevant  to  engineering  design  —  can  very  effectively 
treat  this  larger  class  of  noncoercive  and  nonlinear  equations.  Problems  ad¬ 
dressed  to  date  include  noncoercive  systems  arising  from  frequency-domain 
treatment  of  propagation  phenomena  (the  Helmholtz  equation  [14,15]);  gen- 
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eralized  eigenvalue  problems  [9];  and  the  Burgers  equation  [14,7].  In  addi¬ 
tion,  we  can  now  compute  bounds  not  only  for  the  engineering  outputs,  but 
also  for  the  sensitivity  derivatives  of  these  outputs  with  respect  to  selected 
design  variables  [6]. 

Most  recently  we  have  extended  the  formulation  to  the  full  steady  in¬ 
compressible  Navier-Stokes  equations.  The  theoretical  development  follows 
from  the  Stokes  and  Burgers  treatment,  however  the  difficulties  of  incom¬ 
pressibility  and  nonlinearity  are  compounded,  and  care  must  be  exercised 
in  order  to  obtain  the  desired  (asymptotic)  bounding  property.  We  present 
simple  but  illustrative  results  for  natural  convection  in  an  enclosure  at  low  to 
moderate  Grashof  number.  We  conclude  with  several  remarks  related  to  fur¬ 
ther  extension  of  the  technique,  for  example  to  the  unsteady  Navier-Stokes 
equations. 
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1  Spectral  Eddy- Viscosity  and  Diffusivity 


For  a  recent  review  on  LES,  the  reader  is  referred  to  [l].  Most  of  large- 
eddy  simulation  studies  have  been  developed  in  physical  space.  They  are 
associated  to  a  low-pass  filter  of  width  Ax,  which  is  applied  to  the  flow 
so  as  to  eliminate  fluctuations  in  the  motions  of  smaller  wavelength.  Then,  a 
subgridscale-tensor  appears,  which  is  generally  evaluated  in  terms  of  an  eddy- 
viscosity  function,  such  as  Smagorinsky's  model.  To  me,  the  major  drawback 
of  the  eddy-viscosity  assumption  is  that  it  assumes  a  3pectral  gap  between 
the  resolved  and  subgrid  scales,  assumption  which  i3  never  fulfilled  In  this 
respect,  the  spectral  eddy-viscosity  idea  is  preferable,  provided  one  can  work 
in  Fourier  space,  which  applies  however  only  to  simple  geometries,  as  will  be 
seen. 

We  assume  that  Navier-Stokes  is  written  in  Fourier  space1.  Let  u 
be  the  spatial  Fourier  transforms  of  the  velocity  field.  We  consider  the  cutoff 
wave  number  kc  =  ttAx~1  ,  and  define  a  sharp  low-pass  filter  by  setting  equal 
to  zero  the  amplitudes  at  wave  vectors  k  such  that  |fc|  >  kc-  The  formalism 
of  spectral  eddy  viscosity  is  due  to  Kraichnan  ([3],  see  also  [4]).  One  writes 
the  evolution  equation  for  the  subgrid  modes  k,  and  associate  a  spectral  eddy 
viscosity  to  nonlinear  triadic  interactions  (k,p,  q )  such  that  at  least  p  or  q  is 
larger  than  kc-  In  fact,  this  eddy  viscosity  i3  determined  at  the  level  of  the 
kinetic-energy  spectrum  evolution  equation  given  by  a  two-point  closure  of 
isotropic  turbulence,  die  EDQNM  theory  (see  [5]).  It  is  found  for  3D  isotropic 
turbulence  in  the  case  of  a  Kolmogorov  subgridscale  spectrum: 


vt{k\kc)  =  0.441  CK-V'Xikfkc) 


E{kc) 


kc  \ 


(1) 


where  Ck  is  the  Kolmogorov  constant,  E(kc)  the  kinetic-energy  spectrum 
at  kc,  and  X(k/kc)  a  non-dimensional  function  equal  to  1  up  to  about 
kfkc  =  1/3,  and  sharply  rising  above  (“plateau-peak”  behaviour).  In  the 
same  way,  an  eddy-diffusivity  may  be  introduced  far  a  passive  scalar,  with 
the  corresponding  turbulent  Prandtl  number.  Within  the  EDQNM  theory, 
both  eddy  coefficients  have  a  plateau-peak  behaviour:  it  is  clear  that 


At  this  level,  this  requires  periodicity  in  the  three  spatial  directions.  We  will  see 
below  how  to  get  rid  of  this  assumption. 
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the  plateau  part  corresponds  to  the  usual  eddy-coeffident3  assumption 
when  one  goes  back  to  physical  space,  so  that  the  “peak”  (Kraichnan  called 
it  “cusp" )  part  goes  beyond  the  scale-separation  assumption  inherent  to  the 
classical  eddy-viscosity  and  diffusivity  concepts.  We  have  carried  out  spectral 
LES  of  decaying  isotropic  turbulence[6][7],  where  a  double  filtering  in  Fourier 
space  across  kc  and  k'c  =  kef 2, 

is  performed.  We  find  that  the  subgrid  energy  transfers2  across  k‘c  are 
equal  to  the  subgrid  transfers  across  kc  plus  “resolved”  subgrid  transfers 
across  k!c  corresponding  to  the  low-pass  field  with  respect  to  kc-  This  is,  at 
an  energetic  level,  a  sort  of  Germano’3  identity(2]  in  spectral  space.  These 
calculations  3how  that  the  resolved  eddy  viscosity  does  possess  the  EDQNM 
plateau-peak  shape,  but  not  the  eddy  difiusivity. 

Let  us  3how  now  an  adaptation  of  the  plateau-peak  model  to  kinetic- 
energy  spectra  oc  k~m  for  k  >  kc,  when  the  exponent  m  is  not  necessarily 
equal  to  5/3.  The  spectral  eddy  viscosity  is  now 

ut{k\kc)  =  0.31Ck--^2V3^ 

Sf- X(k/kc)[E(kc)lkc}l/ 2  ,  (2) 

for  m  <  3.  This  expression  is  exact  for  k  «kc  within  the  EDQNM  theory, 
as  shown  in  [7].  We  retain  the  peak  shape  through  X(k/kc)  in  order  to  be 
consistant  with  the  Kolmogorov  spectrum  expression  of  the  eddy  viscosity. 
For  m  >  3,  the  scaling  is  no  more  valid,  and  the  eddy  viscosity  will  be  set 
equal  to  zero.  In  the  spectral-dynamic  model,  the  exponent  m  i3  determined 
through  the  LES  with  the  aid  of  least-3quares  fits  of  the  kinetic-energy  spec¬ 
trum  close  to  the  cutoff,  and  the  turbulent  Prandtl  number  is 

Pr  =  0-18(5  -  m)  (3) 

(see  [5]). 

We  have  applied  the  spectral-dynamic  model  to  the  temporal  mixing  layer 
in  the  case  of  an  initial  3D  white-noise  perturbation,  with  statistical  data 
concerning  velocity,  rms  velocity  fluctuations  and  Reynolds  stresses  in  very 
good  agreement  with  the  experiments  of  unforced  mixing  layers  carried  out 
in  Stanford. 


1.1  Plane  Poiseuille  Flow 

We  present  now  spectral-dynamic  model  results  for  a  periodic  channel.  We 
use  a  mixed  spectral-compact  code,  the  compact  scheme  being  employed  in 
the  transverse  direction,  while  pseudo-spectral  methods  are  used  in  the  longi¬ 
tudinal  and  span  wise  directions  which  are  periodic.  The  channel  has  a  width 
2h,  and  we  define  the  macroscopic  Reynolds  number  by  Re  =  2 hUm/v,  where 
Um  i3  the  bulk  velocity.  The  kinetic-energy  spectrum  allowing  to  determine 

of  momentum,  or  scalar 
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the  eddy-viscosity  13  calculated  at  each  time  3tep  by  averaging  in  planes  par¬ 
allels  to  the  walls,  and  is  thus  a  function  of  ( y,t ).  We  show  now  a  LES  at 
Rt  =  14000  (/i+  =  389).  There  is  a  grid  refinement  close  to  the  wall,  in  order 
to  simulate  accurately  the  viscous  sublayer  (first  point  at  y+  =  1).  We  have 
compared  the  calculation  with  a  DNS  at  h*  =  395  carried  out  by  [8].  Figure 
1  shows  the  mean  velocity  and  the  rms  velocity  components.  The  agreement 
i3  very  good,  which  is  a  severe  challenge  for  the  model.  Notice  that  the  LES 
allows  to  reduce  the  computational  cost  by  a  factor  of  the  order  of  100,  which 
is  huge.  We  stress  also  that  we  have  applied  the  spectral-dynamic  model  to 


Fig.  1.  turbulent  channel  flow,  comparisons  of  the  spectral-dynamic  model  (straight 
lines,  h  =  389)  with  the  DNS  of  Kim  ([8],  symbols,  h  =  395);  a)  mean  velocity, 
b)  nn3  velocity  component3(courtesy  E.  Lamballais). 

the  rotating  channel,  with  good  results  concerning  in  particular  the  linear 
velocity  profile  on  the  anticyclonic  side. 

2  Return  to  Physical  Space 

Now,  let  us  consider  the  EDQNM  eddy  viscosity  (still  scaling  on  \/E(kc)/kc) 
with  no  cusp,  and  adjust  the  constant  by  subgrid-energy  conservation  argu- 
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ment3.  We  assume  that  E{kc,x)  is  a  local  kinetic-energy  spectrum,  calcu¬ 
lated  in  terms  of  the  local  second-order  velocity  structure  function  of  the 
filtered  field 

F2(x,  Ax)  =  <||tt(*,  t)  -  u(x  +  r,  t)||2)||r|1=4s  (4) 

as  if  the  turbulence  is  three-dimensionally  isotropic.  This  yields  for  a  Kol¬ 
mogorov  spectrum 

ufF{x,Ax)^Q.W5C^/2  Ax[F2(x,Ax)}^  .  (5) 

Fz  is  calculated  with  a  local  statistical  average  of  square-velocity  differences 
between  x  and  the  six  closest  points  surrounding  x  on  the  computational 
grid.  In  some  cases,  the  average  may  be  taken  over  four  points  parallel  to  a 
given  plane.  Notice  also  that  if  the  computational  grid  is  not  regular  (but  still 
orthogonal),  interpolations  of  (5)  have  been  proposed  by  [1].  Such  a  structure- 
function  model  (SF)  works  very  well  for  decaying  isotropic  turbulence,  where 
it  yields  a  fairly  good  Kolmogorov  spectrum  ([7]),  better  than  Smagorinsky’s 
model  (with  Cs  —  0.2)  and  the  plateau-peak  models  (simple  or  dynamic). 
However,  it  does  not  work  for  transition  in  a  boundary  layer  at  low  Mach  (or 
incompressible)  where,  like  Smagorinsky,  it  is  too  dissipative  and  prevents 
secondary  instabilities  of  TS  waves  to  develop.  To  overcome  this  difficulty, 
two  improved  versions  of  the  SF  model  have  been  developed.  In  the  selec¬ 
tive  structure-function  model  (SSF),  the  eddy  viscosity  is  set  equal  to  zero  if 
the  angle  between  the  local  vortidty  vector  and  the  average  next-neighbour’s 
vorticity  is  smaller  than  the  value  of  20°,  found  as  the  most  probable  one 
in  equivalent  LES  of  isotropic  turbulence.  Results  concerning  the  coherent 
vortices  in  an  incompressible  badatep  at  Reynolds  5100  will  be  given  at  the 
conference.  We  3how  that  the  Kelvin-Helmholtz  like  vortices  shed  behind 
the  step  undergo  helical  pairing,  and  transform  into  big  staggered  A  vor¬ 
tices  which  impinge  the  bottom  wall  and  are  carried  away  downstream.  We 
determine  the  pressure  spectra  in  terms  of  the  Strouhal  numbers3  at  vari¬ 
ous  locations  downstream  of  the  step.  It  turns  out  that  three  characteristic 
Strouhal  numbers  are  important:  the  harmonic  0.23,  the  subharmonic,  and 
the  Sapping  frequency  of  the  recirculating  bubble  0.07. 

In  the  filtered  structure-function  model  (FSF)[10],  the  filtered  field  U{  is 
submitted  to  a  high-pass  filter  in  order  to  get  rid  of  low-frequency  oscillations 
which  affect  E„  ( kc )  in  the  SF  model.  The  high-pass  filter  is  a  Laplacian  dis¬ 
cretized  by  second-order  centered  finite  differences  and  iterated  three  times. 
It  is  found 

ufSF(x,  Ax)  =  0.0014  C~3/2  Ax  [/2(®,  Ax)]1'2  .  (6) 


non-dimeasionalized  by  the  incoming  velocity  and  the  step  height 
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3  Compressible  Boundary  Layer 

We  present  here  some  recent  results  obtained  in  Grenoble  concerning  LES 
of  transition  to  turbulence  in  a  boundary  layer  above  an  adiabatic  flat  plate 
in  an  ideal  gas  at  =  0.3.  We  have  developed  a  formalism  using  density- 
weighted  Favre  averages  [9].  One  introduces  a  “macro-temperature”  and  a 
“macro-pressure”  which  are  related  by  the  equation  of  state 

xa  ~  pRd  ...  (7) 

provided  the  condition 

5  7^3gs  <  1 

is  satisfied.  For  the  air,  it  improves  the  condition  lAffgs  <  1  proposed  usu¬ 
ally.  With  the  classical  eddy  assumptions  done  for  the  momentum  and  heat 
subgrid  stresses,  it  permits  to  reduce  the  problem  to  the  resolution  of  com¬ 
pressible  Navier-Stokes  equations,  where  the  molecular-diffusion  coefficients 
for  momentum  and  heat  are  supplemented  by  incompressible  eddy  counter¬ 
parts.  Here,  the  resolution  at  the  wall  is  y4"  =  1  or  2,  and  the  Mach  number 


-  cf  y’l 

- cf  y*2 

cf  coustots 
•  cf  tarenbkitt 


Fig.  2.  LES  of  a  spatial  boundary  layer  at  Mach  0.3:  friction  coefficients  against 
downstream  distance,  compared  with  theoretical  predictions  of  Cousteix  and  Baren- 
blatt  (courtesy  E.  Briand). 


0.3.  Upstream  conditions  (harmonic  K-mode  or  3ubharmonic  H-mode)  are 
obtained  with  the  aid  of  nonlinear  PSE  calculations(ll]  using  the  ONERA 
code[12|.  To  this  upstream  state  (with  still  R$  —  1000),  one  superposes  a 
3D  white-noise  of  amplitude  0.2  the  amplitude  of  the  PSE  perturbation.  The 
simulations  involve  up  to  five  millions  of  grid  points.  Figure  2  shows  for  the 
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K-caae  the  downstream  evolution  of  the  friction  coefficient  at  the  wall,  with 
comparison  against  the  theoretical  predictions  of  Van  Driest  and  Barenblatt. 
One  sees  a  good  agreement  of  the  LES  with  these  predictions,  with  an  im¬ 
provement  with  a  resolution  of  y+  =  1.  It  is  even  better  in  the  H-case. 


3.1  Boundary  Layer  upon  a  Groove 

We  give  here  results  concerning  a  turbulent  boundary  layer  on  a  flat  plate  at 
Rs  «  1000  passing  over  a  square  2D  groove.  The  Mach  number  is  0.5,  and 
the  flow  upstream  is  calculated  through  an  iterative  procedure  developed  by 
[13].  The  cavity  has  a  depth  equal  to  the  upstream  boundary-layer  thickness 

S.  We  have  looked  at  the  vorticity  modulus  map  conditioned  by 

the  “Q*criterion."[14],  which  selects  regions  where  the  second  invariant  of 
the  velocity-gradient  tensor  is  positive.  The  calculation  shows  how  boundary- 
layer  hairpins  shed  downstream  of  the  backward  step  are  lifted  above  the 
recirculation  zone  within  the  cavity,  and  impinge  the  upward  step  to  form 
both  smaller-3cale  turbulence  (-in  a  sort  of  ultraviolet  kinetic-energy  cascade) 
and  big  inclined  arch-like  vortices  which  travel  downstream  of  the  cavity.  One 
can  check  also  that  the  characteristic  3p  an  wise  wavelengths  of  the  velocity 
streaks  decrease  after  the  cavity. 
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Summary 

Spectral  methods  have  been  used  with  great  success  in  direct  numerical  simulations  of  turbulent  flows 
m  simple  computational  domains.  In  this  talk  we  will  present  the  next  generation  of  spectral  methods 
on  unstructured  and  polymorphic  domains  suitable  for  addressing  the  geometric  complexity  of  industrial 
applications.  These  new  methods  are  hierarchical  and  are  well  suited  for  computational  steering  and  dy¬ 
namic  direct  numerical  simulation  (dDNS).  We  will  present  several  dDNS  examples  of  both  internal  and 

open  transitional  and  turbulent  flows  to  demonstrate  this  new  methodology  and  compare  it  with  low-order 
methods. 

hi  particular,  we  will  present  a  detailed  study  of  flow  past  a  flexible  cylinder  subject  to  vortex  induced 
vibrations  (VIV).  The  question  of  what  flow  pattern  prevails  in  the  wake  of  a  long  flexible  structure  sub¬ 
ject  to  uniform  or  shear  flow  is  of  fundamental  importance  as  it  dictates  the  magnitude  of  the  forces  on 
the  structure..  What  conditions  produce  a  stable  standing  wave  response  or  a  travelling  wave  response? 
We  have  obtained  simulation  results  and  models  from  this  ongoing  investigation  using  three-dimensional 
simulations  of  flow  past  a  flexible  cable.  Different  initial  conditions  and  different  spans  of  the  cable  were 
used  to  classify  the  responses.  Travelling  waves  appear  in  long  cables  and  produce  oblique  shedding  with 
a  wake  similar  to  the  two-dimensional  von  Karman  street  across  the  entire  span.  In  contrast,  standing 
waves  appear  in  shorter  cables  and  produce  parallel  shedding  with  a  wake  structure  varying  continuously 
across  the  span,  from  a  von  Karman  street  at  the  anti-nodes  to  anti-symmetric  shedding  at  the  nodes. 
Chevron-like  shedding  is  produced  in  the  presence  of  shear  along  the  cable. 

Another  important  question  for  VIV  is  the  wake  transition  process.  Transition  to  turbulence  in  the 
wake  of  a  fixed  cylinder  occurs  at  Reynolds  number  between  250  and  400,  as  has  been  established  by 
experimental  results  and  our  previous  simulation  studies.  However,  the  transition  process  changes  funda¬ 
mentally  if  the  cylinder  is  flexible,  and  vibrates  freely.  In  this  numerical  study,  we  considered  flow  past 
exible  beams  and  cables  undergoing  free  oscillations  subject  to  near  lockin  excitation.  We  investigated 
i  erent  vibrating  conditions,  corresponding  to  varying  the  bending  stiffness  (beams),  tension  (cables)  and 
the  mass  ratio  parameter,  in  order  to  determine  the  new  transition  mechanisms.  In  general,  cables  tend 
to  promote  wake  transition  whereas  beams  tend  to  delay  transition  compared  to  the  fixed  cylinder  behavior. 
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Abstract.  A  Large  Eddy  Simulation  of  the  flow  over  a  matrix  of  surface  mounted 
cubes  is  presented.  Computations  employing  different  sub-grid  scales  models  and 
different  domain  sizes  are  performed.  Mean  velocity  and  Reynolds  stresses  profiles 
are  evaluated  and  compared  with  experiment  and  with  previous  single-cube  LES 
computations.  The  influence  of  the  periodic  boundary  conditions  is  discussed  by 
doubling  the  size  of  the  sub-channel  domain  considered  for  the  simulation. 


1  Introduction 


The  turbulent  flow  around  bluff  bodies  is  characterized  by  complex  inter¬ 
actions  between  different  phenomena  such  as  boundary  layers,  shear  layers, 
separation  and  reattachment.  Such  flows  are  relevant  for  many  environmen¬ 
tal  and  industrial  applications.  Although  statistical  turbulence  models  have 
been  applied  successfully  in  many  practical  computations,  the  Large  Eddy 
Simulation  technique  shown  proved  its  better  potential  in  calculating  flows 
characterized  by  large  unsteady  vortical  structures.  In  particular  when  the 
mixing  of  a  scalar  such  as  temperature  is  of  concern,  the  unsteady  flow  has 
to  be  determined  with  high  accuracy.  It  then  allows  to  detect  for  example 
local  overheating  of  the  surface  and  may  furnish  relevant  information  for  the 
design  of  cooling  devices,  burners,  etc.  In  the  present  study,  a  Large  Eddy 
Simulation  of  a  flow  over  am  array  of  surface-mounted  cubical  obstacles  is 
considered.  This  flow  is  an  example  of  a  bluff  body  flow  where  large  scale 
structures  dominate  the  turbulent  transport,  and  for  which  detailed  experi¬ 
mental  results  [1]  are  available. 

2  Numerical  Method 

The  Large  Eddy  Simulation  code  LESOCC  develpped  at  the  Institute  for 
Hydromechanics  [2,3]  is  based  on  a  Finite- Volume  method  for  solving  the 
incompressible  Navier-Stokes  equations  on  general  body-fitted,  curvilinear 
grids  (LESOCC  =  Large  Eddy  Simulation  On  Curvilinear  Coordinates).  A 
non-staggered,  cell-centered  grid  arrangement  is  used.  Both  convective  and 
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viscous  fluxes  are  approximated  by  central  differences  of  second  order  accu¬ 
racy.  The  temporal  discretization  consists  of  an  explicit  low-storage  Runge- 
Kutta  method  (second  order  in  time).  The  conservation  of  mass  is  achieved 
by  a  standard  pressure  correction  algorithm  (SIMPLE).  The  Poisson  equa¬ 
tion  for  the  pressure  correction  variable  is  solved  by  the  strongly  implicit 
procedure  of  Stone  [4],  and  in  order  to  avoid  decoupling  of  pressure  and  ve¬ 
locity  on  the  non-staggered  grid  the  momentum  interpolation  proposed  by 
Rhie  and  Chow  [5]  is  applied.  The  code  is  highly  vectorized  and  has  been 
validated  extensively,  also  for  the  case  of  a  single  surface  mounted  cube  at 
two  Reynolds  numbers  [3,6].  Different  subgrid-scale  models  are  implemented: 
the  Smagorinsky  model  with  van  Driest  damping  near  solid  walls,  and  the  dy¬ 
namic  model  with  the  least-squares  approach  of  Lilly.  Different  wall  function 
models  are  also  implemented  but  are  not  employed  in  the  present  compu¬ 
tations  because  the  Reynolds  number  is  fairly  low.  Recently  the  code  has 
been  extended  to  block-structured  grids  for  parallel  execution  on  distributed 
memory  machines. 

3  Details  of  the  test-case  calculations 

An  experimental  situation  is  simulated  in  which  a  matrix  of  cubes  of  height 
H  is  placed  on  one  of  the  walls  of  a  two-dimensional  channel  of  height  3.4#  as 
sketched  in  Figure  I.  The  Reynolds  number  based  on  the  height  of  the  channel 
and  the  bulk  velocity  is  Re  =  13000.  Due  to  the  large  number  of  equispaced 
cubes  with  distance  S  =  4 #,  the  mean  flow  is  periodic  with  a  period  equal  to 
4if  in  the  x-  and  z-direction  [1],  A  sub-channel  of  size  AH  x  3.4#  x  4#  (in  x-, 
y-  and  z-direction  respectively)  is  considered  in  the  simulation,  using  periodic 
boundary  conditions  in  the  streamwise  and  in  the  spanwise  direction. 


Fig.  1.  Sketch  of  the  matrix  configuration  of  cubes  on  the  channel  wall  (£rom[l]) 
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To  assess  the  influence  of  the  periodic  boundary  conditions  in  the  stream- 
wise  direction,  a  computation  was  performed  by  doubling  the  size  of  this 
sub-channel  domain.  All  the  computations  were  performed  on  a  stretched 
grid  with  100  x  100  x  100  and  200  x  100  x  100  grid  points  for  the  x-,  y-  and 
z-  direction,  for  the  standard  and  doubled  domain  respectively.  For  the  par¬ 
allel  execution,  the  standard  domain  was  divided  into  4  blocks,  with  2  blocks 
in  both  x-  and  z-direction,  and  the  second  one  was  divided  into  8  blocks. 
In  this  paper  four  different  computations  are  presented.  Three  computations 
were  performed  with  the  standard  domain,  employing  the  Smagorinsky  model 
(SRUN4),  the  dynamic  model  (DRUN4),  and  no  sub-grid  model  (NRUN4). 
One  computation  (DRUN8)  was  performed  employing  the  doubled  domain 
size  and  the  dynamic  model. 


Table  1.  Simulations 


RUN 

Comp.  Domain 

Grid  size 

Model 

NRUN4 

4iT  x  3.41?  x  4J? 

No  SGS  Model 

SRUN4 

4 if  x  3.4J2-  x  4 if 

Smagorinsky 

DRUN4 

4 3  x  3.4if  x  4 if 

Dynamic 

DRUN8 

8 if  x  Z.iff  x  iff 

Dynamic 

4  Results  of  the  simulation 

4.1  Description  of  the  general  flow  structure 


From  visualization  studies  and  detailed  measurements,  Meinders  et  al.  [1]  de¬ 
vised  a  sketch  of  the  three-dimensional  flow  pattern  around  the  cube  for  one 
spatial  period  in  the  matrix.  In  this  part  we  will  discuss  the  comparison  be¬ 
tween  these  experimental  results,  the  numerical  simulation  and  the  previous 
experiment  [3,6]  and  simulations  of  the  flow  around  a  single  cube. 

The  flow  is  characterized  by  the  presence  of  distinct  vortex  structures  in 
the  vicinaty  of  the  obstacles.  As  shown  in  Figure  2,  the  mean  flow  separates  in 
front  of  the  cube,  generating  a  primary  and  a  secondary  vortex.  The  primary 
vortex  is  bent  as  a  horseshoe  vortex  around  the  cube  into  the  wake.  The  flow 
separates  also  at  the  front  corners  of  the  cube  on  the  roof  and  side  walls.  A 
large  separation  region  develops  behind  the  cube,  with  the  development  of 
an  arch  vortex  interacting  with  the  horseshoe  vortex.  This  global  structure 
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is  close  to  the  one  observed  for  a  single  cube,  however  due  to  the  confine¬ 
ment  of  the  flow  in  the  spanwise  direction,  and  due  to  the  interaction  in  the 
streamwise  periodic  direction  between  neighboring  horseshoe  type  vortices, 
the  reattachment  on  the  roof  and  on  the  side  walls  occurs  at  a  different  posi¬ 
tion.  In  the  matrix  case,  the  reattachment  length  on  the  roof  and  on  the  side 
walls  is  approximately  0.3ff ,  whereas  the  single-cube  calculations  showed  no 
reattachment  on  the  roof.  Nevertheless  the  reattachment  length  behind  the 
cube  is  very  similar  in  both  the  single  and  matrix  cube  case  (approximately 

1.5  cube  heights  from  the  leeward  face). 


Fig.  2.  Streamlines  of  the  time-averaged  flow:  in  the  symmetry  plane  and  on  the 
floor  of  the  channel  (data  from  SRUN4  duplicated  in  the  streamwise  direction) 
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4.2  Vortex  shedding 

As  for  the  single-cube  configuration,  the  experiment  showed  a  predominant 
fluctuation  frequency  sideways  behind  the  cube,  corresponding  to  vortex 
shedding  as  the  flow  passes  the  side  walls.  These  coherent  structures  were  de¬ 
tected  from  power  density  spectra  of  the  time  series  of  the  spanwise  velocity 
component.  The  dominant  characteristic  frequency  derived  from  the  location 
in  the  spectrum  which  corresponded  to  the  maximum  energy  around  27 Hz 
leading  to  a  Strouhal  number  St  =  SL  -  0.109.  The  same  procedure  applied 
to  the  Large  Eddy  Simulation  gives  a  Strouhal  number  St  =  0.12  which  is 
close  to  this  value. 

4.3  First  and  second  order  moments 


Detailed  comparison  between  'experimental  and  numerical  mean  velocities 
and  Reynolds  Stresses  at  different  locations  were  the  subject  of  a  recent 
Workshop  [7].  A  set  of  these  profiles  is  given  in  Figure  3.  While  streamwise 
fluctuations  u/2  become  largest  in  the  shear  layer  and  in  front  of  the  cube 
(r/h  =  3.7),  the  3panwise  fluctuations  become  maximum  close  to  the  reat¬ 
tachment  point  ( zjh  s  2.3).  The  overall  agreement  between  the  experiment 
and  the  simulation  is  fairly  good,  despite  a  slight  discrepancy  in  the  wake 
for  the  normal  Reynolds  Stress  w2  which  has  been  also  observed  in  similar 
computations  [8].  The  different  sub-grid  scale  models  give  almost  the  same 
results,  with  slightly  better  results  for  SRUN4.  They  do  not  significantly 
improve  the  results  obtained  with  NRUN4  (without  sub-grid  model).  The 
simulation  conducted  with  two  cubes  included  in  the  computational  domain 
does  not  show  significant  difference  with  the  other  simulations.  Hence  inter¬ 
actions  between  consecutive  horseshoe  vortices  have  only  very  little  influence 
on  first  and  second  order  moments. 

5  Computing  times 


All  computations  were  performed  on  the  Vector  Parallel  mainframe  Fujitsu 
VPP  300  installed  at  the  Computer  Center  of  the  University  of  Karlsruhe. 
It  has  16  processors  each  having  a  peak  performance  of  2.2  Gfiops  and  a 
core  memory  of  2  Gbytes.  Typical  CPU  times  per  time  step  and  grid  point, 
and  speed-up  during  parallel  computations  are  given  Table  1.  The  result  of 
a  computation  using  only  one  block  (DRUN1)  is  also  included  in  the  table. 
The  presented  results  were  obtained  with  a  time  step  &t  —  2.2  x  10“3.  15000 
steps  were  computed  for  the  startup  phase  and  statistics  were  acumulated 
during  90000  time  steps.  The  whole  runs  took  about  50  hours  of  CPU  time 
per  processors. 
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Fig.  3.  Mean  velocities  and  Reynolds  stress  profiles  at  different  positions  in  the 
streamwise  direction:  x/H  -  0.3;  x/H  =  1.3;  x/H  =  1.7;  z/H  =  2.3;  x/H  =  2.7; 
x/H  =  3.7 
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Table  2.  Computation  times 


RUN 

Nb  of  points 

Nb  of  PE 

CPU/tim*  •t«p/fri4  poi«n 

Speed-up 

DRUNl 

1.10* 

1 

9.  10"* 

— 

DRUN4 

1.10* 

4 

2.39  10"* 

3.7 

DRUN8 

2.10* 

8 

1.21  10"* 

7.4 
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Traditional  Large-Eddy  Simulation  (LES)  provides  good  results  on  standard  benchmark  turbulent  flows, 
but  engineering  applications  of  more  complex  geometries  are  restricted.  This  is  mainly  due  to  the  phys¬ 
ical  hypothesis  of  homogeneity,  isotropy,  or  local  equilibrium  of  the  subgrid  scale  models,  which  are 
hazardous  in  more  realistic  configurations.  In  an  attempt  to  overcome  these  problems,  the  recent  notion 
of  Very  Large-Eddy  Simulation  (VLES)  has  been  developed.  LES  relies  on  an  equation  filtering  theory, 
actually  governed  by  the  modeling  terms  whose  influence  on  the  flow  dynamics  is  not  well  described.  It 
appears  essential  to  qualify  these  built-in  filters-of  LES  models  before  performing  computations  on  coarse 
mesh/high  Reynolds  number.  This  work  is  an  assessment  of  previous  theoretical  studies  on  LES  filters 
and  kinetic  energy  spectrums. 

In  the  present  study,  computations  have  been  carried  out  on  homogeneous  isotropic  turbulent  ,  and  freely 
decaying  flow.  The  initial  random  field  has  a  gaussian  kinetic  energy  spectrum  centered  on  a  wavenumber 
ko  =  2.  The  discretization  is  (66)^  points  in  a  periodic  cube  of  side  2 <r. 

The  PEGASE  code,  developed  at  ONERA,  solves  the  full  unsteady  compressible  Navier-Stokes  equations 
using  conservative  variables  (p,  pui,  pEj)  with  a  finite  difference  scheme  on  structured  uniform  grids. 
Convection  terms  are  written  in  skew-symmetric  form  and  along  with  diffusion  terms,  are  discretized 
using  a  fourth-order  centered  scheme.  Time  integration  is  performed  with  a  third-order  low-storage 
Runge-Kutta  scheme.  The  well-known  Smagorinsky  model  [6],  and  the  Schumann  model  [5]  (a  subfilter 
viscosity  depending  on  the  subfilter  kinetic  energy  plus  a  transport  equation  based  on  Prandtl-Kolmogorov 
equation)  have  been  used.  Mach  number  is  chosen  0.1,  so  compressiblity  has  no  significant  effect,  and 
molecular  viscosity  is  null  (Re  =  oo). 

A  similarity  theory  of  Smagorinsky-type  LES,  has  been  proposed  by  Muschinskv  [4].  With  A  bein<^  the 
filter  length  and  A  the  mesh  width  we  define  the  dimensionless  ratio  r  =  A/A  and  the  filter  cut-off  wave 

number  kc  =  it/ A.  Then  the  turbulent  kinetic  energy  (TKE)  spectrum  generated  by  a  ‘Smagorinsky 
fluid’  may  be  written: 

E(k )  =  *o css(^2/3k-^fLES(k/kc,  r) 

where  e  is  the  dissipation,  k  the  wave  number  and  fLBS  is  a  damping  function.  The  filter  is  mathemat- 
ica  Iy  defined  in  the  physical  space  as  a  convolution  integral  similar  to  a  multiplication  in  the  spectral 
space.  Thus,  in  a  filtering  viewpoint,  we  can  identify  fLBS  as  the  square  transfer  function  G2. 

The  subfilter  viscosity  of  the  Smagorinsky  model  is  uLBS  =  (C5A)2|S|  where  |S|  is  the  modulus  of  the 
deformation  tensor,  the  Smagorinsky  constant  being  Cs  =  0.18.  Whatever  the  form  of  E(k)  it  is  shown 
that  the  dissipation  length  for  a  ‘Smagorinsky  fluid’  is  r,LBS  =  (^„/e)l/4  =  C5A.  A  series  of  simulations 
have  been  carried  out  using  values  of  r  ranged  from  1,  the  usual  value  in  LES,  to  4.  For  a  given  mesh 
size  this  corresponds  to  a  change  of  the  subfilter  scale  (fig.  1).  It  is  shown  that  (fig.  2  left): 

•  We  find  that  for  large  ratio  r  the  damping  function  best  match  the  Heisenberg  [1]  function: 

fff(z)  =  (H-xV/3,  x  =  k/kc 

and  that  the  length  A  controls  the  filter  cut-off  wavenumber.  But  for  r  =  1,  the  filter  is  a  sharp 
cut-off  filter  and  /cbs(x)  =  H(l  —  x)  (H  is  the  Heaviside  function). 
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Simulations  on  a  (126)3  grid,  with  different  initial  conditions,  or  with  a  statistically  steady  forced  turbu¬ 
lence  provide  similar  results. 

Simulations  have  also  been  performed  for  the  Schumann  model,  an  other  eddy-viscosity  model: 

VLE5  =  CAy/  C  —  0.09 

,  2  2  3/2 

Jg£  =  tles:S  -  C^-  +  C2V.(A^Vq2Sr)  +^qh  Ci  =  1,  C2  =  0.1 

where  q2F  is  the  subfilter  kinetic  energy,  and  tlbs  is  the  Reynolds  stress  tensor.  In  this  case  the  damping 
function  has  an  exponential  form  (fig.l  right),  and  the  slope  of  /n(£(£))  depends  linearly  on  the  ratio  r. 
In  terms  of  dimensionless  parameters  the  damping  function  is 

fiEs  =  ea(r)+0(p)*f  x  =  k/kCl  /3(r)  =  — -4.5 

r 

A  similar  exponential  form  of  the  energy  spectrum  in  the  near-dissipation  range  has  been  found  by  pre¬ 
vious  direct  numerical  simulations  [2]. 

For  the  Smagorinsky  model  the  Kolmogorov  ‘constant’  Ko(r)  is  determined  by  ploting  E(k)e~2^3k5^ffss 
(fig.2  left).  From  these  simulations  we  conclude  that: 

•  Kqles  depends  on  r  as  predicted  [4],  considering  the  given  resolution,  and  the  related  turbulent 
Reynolds  number  Rclbs  —  (N/ A)4/3  (fig.2  right). 
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Figure  1:  Left:  TKE  spectrum  of  Smagorinsky-LES  for  ratio  r  =  1  (bold  line),  1.5,  2,  2.5,  3,  3.5, 
and  4.  Right:  Compensated  energy  spectrum  £(A:)e“2/3A3/3  for  the  Smagorinsky  (dashed  line)  and  the 
Schumann  (bold  line)  models  for  r  =  1,2, 2.5, 3, 3.5, 4. 


Figure  2:  Left:  Spectrum  of  £’(fc)e-2/3Ars/3/-i  for  r  =  2  (solid  line),  r  =  3  (dashes),  r  =  4  (dots),  and 
compensated  energy  for  r  =  1  (symbols).  Right:  Kolmogorov  constant  vs.  ratio  r  computed  from  freely 
decaying  homogeneous  and  isotropic  turbulence  of  a  ’Smagorinsky  fluid’. 
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ABSTRACT 

What  are  the  key  issues  for  successful  applications  of  the  large  eddy  simulation  (LES) 
technique  to  practically  relevant  flow  problems  ? 

•  First,  technical  applications  are  in  general  concerned  with  complex  geometries  which 
can  either  be  tackled  by  unstructured  grids  or  by  block-structured  curvilinear  body- 
fitted  grids.  All  numerical  methods  based  on  Cartesian  co-ordinates  or  with  specific 
restrictions  such  as  periodicity  in  a  certain  direction  (e.g.  FFT)  fail  in  this  context. 

•  Second,  industrial  applications  are  in  most  cases  high  Reynolds  number  flows.  As 
known  from  the  considerations  of  Kolmogorov,  the  ratio  of  the  largest  to  the  small¬ 
est  length  scales  ( L/lk )  strongly  increases  with  the  Reynolds  number  leading  to  a 
broader  spectrum  of  turbulent  eddies.  Consequently,  LES  of  high  Reynolds  number 
flows  make  special  demands  with  respect  to  the  applied  numerical  methods  and  the 
subgrid  scale  (SGS)  models.  Assuming  that  the  size  of  the  computational  grid  can¬ 
not  be  enlarged  according  to  Re9/4,  the  SGS  model  has  to  take  a  wider  spectrum 
of  turbulent  vortices  into  account.  Another  topic  of  rising  complexity  for  higher 
Reynolds  numbers  is  the  formulation  of  appropriate  boundary  conditions  especially 
at  solid  walls. 

•  Third,  because  it  is  well-known  that  LES  is  not  a  cheap  tool  for  computing  turbulent 
flows,  its  application  should  be  restricted  to  flow  problems  for  which  an  appropri¬ 
ate  description  by  Reynolds-averaged  Navier-Stokes  equations  (RANS)  combined 
with  statistical  turbulence  models  is  difficult  to  achieve.  An  illustrative  example 
represents  the  flow  past  bluff  bodies  which  is  in  general  very  complicated,  including 
complex  phenomena  such  as  separation,  reattachment  or  vortex  shedding. 

The  long-term  objective  of  the  work  reported  here  is  to  develop  a  LES  technique  which  is 
able  to  simulate  high  Reynolds  number  flows  of  practical  relevance,  especially  bluff  body 
flows.  With  respect  to  the  considerations  above,  the  developed  code  ( CSSQCC  =  £arge 
£ddy  Simulation  On  Curvilinear  Coordinates)  is  based  on  a  3-D  finite-volume  method  for 
arbitrary  non-orthogonal  and  non-staggered  grids.  Details  about  the  temporal/spatial 
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discretization  and  the  implemented  SGS  models  (Smagorinsky  and  dynamic  model)  can 
be  found  in  [1,3,4, 6j.  Recently  C6SOCC  has  been  extended  by  a  multi-block  structure 
strongly  improving  the  possibility  to  resolve  complex  geometries.  Furthermore  the  multi¬ 
block  implementation  was  also  the  basis  for  parallelization  by  domain  decomposition  and 
message  passing  (MPI).  C6SOCC  is  highly  vectorized  (vectorization  rate  >  99.8%)  allow¬ 
ing  to  perform  efficient  computations  on  vector-parallel  machines  such  as  NEC  SX-4  or 
Fujitsu  VPP  300/700.  The  simulations  presented  in  this  abstract  have  been  performed 
mainly  on  a  VPP  300  using  4  processors. 

Before  CSSOCC  was  applied  to  the  flow  problem  of  the  present  investigation,  it  was  ex¬ 
tensively  tested  and  verified  by  LES  computations  of  a  variety  of  different  test  cases.  This 
included  plane  channel  flow,  flows  through  a  straight  square  duct  and  a  180°  bend  [1], 
flow  around  a  surface-mounted  cubic  obstacle  placed  in  a  plane  channel  [2, 4, 8]  and  flow 
past  a  long  square  cylinder  [3,8].  In  the  course  of  this  validation  process  detailed  inves¬ 
tigations  on  different  numerical  and  modeling  aspects  influencing  the  quality  of  the  LES 
solution  have  been  performed  for  the  flow  past  a  long  circular  cylinder  at  a  low,  subcritical 
Reynolds  number  of  Re  —  3900  [5,6].  In  contrast  to  its  square  counterpart,  this  flow  con¬ 
figuration  requires  curvilinear  body-fitted  grids.  Vortex  shedding  past  a  circular  cylinder 
is  also  physically  more  challenging  than  the  square  cylinder  test  case  because  the  separa¬ 
tion  point  on  the  surface  of  the  cylinder  is  not  fixed  by  the  geometry.  These  characteristics 
of  the  circular  cylinder  flow  combined  with  the  enormous  number  of  experimental  studies 
on  the  vortex  dynamics  of  cylinder  wakes  (see,  e.g.,  review  by  Williamson  [9])  reveals 
this  flow  to  be  an  excellent  test  case  for  the  intended  investigations  on  LES  of  bluff  body 
aerodynamics.  Figure  1  shows  the  instantaneous  von  Kannan  vortex  street  past  the  cylin¬ 
der  visualized  by  streaklines  and  the  time-averaged  streamlines.  For  this  low  Reynolds 
number  case  {Re  =  3900)  good  agreement  was  found  with  experimental  measurements 
available.  For  details  see  [6]. 

Concerning  the  considerations  mentioned  above,  these  investigation  were  now  extended 
to  the  more  challenging  test  case  of  circular  cylinder  flow  at  Re  —  1.4*  10s  which  has  been 
studied  experimentally  by  Cantwell  and  Coles  [7].  At  this  higher  Reynolds  number  the 
flow  is  still  subcritical  so  that  transition  takes  place  in  the  free  shear  layers.  Two  different 
curvilinear,  O-type  grids  have  been  generated.  The  first  consists  of  165  x  165  control 
volumes  (CV)  in  the  cross-sectional  plane  and  64  computational  cells  in  the  spanwise 
direction.  The  second  is  refined  in  the  cross-sectional  plane  using  325  x  325  CV  while 
the  spanwise  resolution  is  kept  constant  (total  of  6.76  Mill.  CV).  The  points  are  clustered 
in  the  vicinity  of  the  cylinder  in  order  to  resolve  the  extremely  thin  boundary  layer  at 
the  cylinder  and  to  apply  no-slip  wall  boundary  conditions.  Additionally,  clustering  of 
grid  points  is  performed  for  the  wake  region.  The  entire  integration  domain  has  a  radial 
extension  of  1 5D  in  the  cross-section  (D  =  cylinder  diameter).  Three  different  spanwise 
extensions  (ttD,  2D,  D )  have  been  investigated.  In  this  direction  of  the  cylinder  periodic¬ 
ity  of  the  flow  is  assumed.  At  the  inflow  plane  constant  velocity  is  imposed.  A  convective 
boundary  condition  is  used  at  the  outflow  boundary.  Statistics  are  compiled  over  several 
vortex  shedding  cycles  and  in  the  spanwise  direction. 

Figure  2  shows  first  qualitative  results  of  the  predicted  flow  past  a  circular  cylinder  at 
Re  —  1.4*10°.  In  comparison  with  the  low  Reynolds  number  case  the  recirculation  region 
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behind  the  cylinder  (see  Fig.  2  (a):  time— averaged  streamlines)  is  much  smaller  which  is 
in  good  agreement  with  experiments  [7].  Figs.  2  (b)  and  (c)  show  contours  of  the  total 
resolved  Reynolds  stresses  u'u'  and  u'v' .  In  the  paper  these  quantities  will  be  compared 
with  the  experimental  values  [7].  Finally,  the  time  histories  of  the  drag  and  lift  coefficient 
(Cd  and  Ci)  are  plotted  in  Fig.  2  (d).  Apart  from  cyclic  oscillations  due  to  the  vortex 
shedding  phenomenon  and  high-frequency  turbulent  fluctuations,  the  signals  clearly  have 
a  low-frequency  component  which  modulates  the  time  histories  of  Cd  and  Ct.  However, 
in  comparison  with  the  low  Reynolds  number  case  (not  shown  here)  the  modulation  of 
the  signals  is  much  weaker  for  Re  =  1.4  *  10s. 

A  more  detailed  description  of  all  aspects  which  have  been  investigated  and  which  play  a 
dominant  role  for  the  quality  of  LES  predictions  will  be  presented  and  discussed  in  the  full 
paper.  This  will  include  a  quantitative  comparison  of  mean  values  as  well  as  turbulence 
quantities.  Experimental  data  available  [7]  will  be  used  for  verifying  the  computed  results 
and  supporting  the  investigations. 
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(a)  Streaklines. 


(b)  Time-averaged  streamlines. 


Figure  1:  Large  eddy  simulation  of  the  flow  past  a  circular  cylinder  at  Re  =  3900. 


(a)  Time-averaged  streamlines. 


(b)  Contours  of  vJu‘ 


Figure  2:  Large  eddy  simulation  of  the  flow  past  a  circular  cylinder  at  Re  —  1.4  *  10°. 
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Numerical  simulation  of  turbulence  has  its  roots  in  the  atmospheric  sci¬ 
ences  and  in  particular  in  Numerical  Weather  Forecasting  (NWP).  Although 
the  first  ideas  on  NWP  were  formulated  by  Richardson  in  1922,  serious 
progress  in  NWP  was  only  possible  after  the  advent  of  electronic  computers 
just  after  the  second  world  war.  As  a  spin  off  from  NWP  the  first  numerical 
simulations  of  turbulent  convection  were  performed  first  in  two  dimensions 
by  Lilly  and  somewhat  later  in  three  dimension  by  Deardorff.  Especially, 
the  work  of  Deardorff  has  become  widely  known  and  it  has  inspired  the 
use  of  numerical  simulation  in  other  areas  of  turbulence  research.  Many  of 
the  techniques  that  Deardorff  introduced  and  the  results  that  he  obtained, 
are  still  relevant  today.  As  example  we  may  mention  here  the  Smagorinsky 
model  which  was  originally  developed  as  a  model  to  parameterize  the  small 
scale  motions  in  large-scale  meteorological  models.  With  respect  to  the  at¬ 
mospheric  boundary  layer  Deardorff  performed  in  the  early  1970’s  the  first 
LES  of  the  so-called  convective  boundary  layer.  This  is  the  turbulent  layer 
which  develops  over  a  land  surface  and  in  which  turbulence  is  generated  as 
a  result  of  natural  convection  due  to  solar  heating.  Since  then  LES  com¬ 
putations  of  the  convective  boundary  layer  have  been  repeated  many  times 
and  one  could  say  that  our  general  knowledge  of  this  type  of  boundary  layer 
is  mainly  due  to  LES. 

As  mentioned  above,  the  use  of  LES  has  been  in  particular  successful 
for  the  convective  boundary  layer.  The  results  of  LES  have  led  here  to 
the  formulation  of  so-called  mixed-layer  scaling.  In  this  scaling  approach 
the  variables  in  the  convective  boundary  when  scaled  in  terms  of  a  limited 
number  of  characteristic  parameters,  collapse  to  universally  valid  scaling 
relationships.  After  mixed-layer  scaling  was  derived  based  on  LES  data,  it 
was  also  confirmed  by  laboratory  and  atmospheric  experiments.  A  second 
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useful  application  of  LES  in  the  atmosphere  is  the  neutral  boundary  layer. 
In  this  type  of  boundary  layer  temperature  effects  are  absent.  This  is  an 
exception  in  the  real  atmosphere  and  therefore  it  is  difficult  to  observe  a 
neutral  boundary  layer.  In  this  case  LES  is  the  only  alternative. 

Experiments  in  the  atmospheric  boundary  layer  are  very  expensive  and 
also  very  time  consuming.  Moreover,  a  completely  reproducible  and  con¬ 
trolled  experiment  in  the  atmosphere  is  impossible  by  definition  because 
of  the  natural  variability  of  the  atmosphere.  Therefore  LES  seems  an  at¬ 
tractive  alternative.  Considering  the  successes  of  LES  in  the  atmosphere 
as  mentioned  above,  it  seems  tempting  the  use  LES  as  an  alternative  for 
experiments.  However,  despite  its  successes  LES  also  knows  some  failures. 
Examples  are  the  behaviour  of  the  horizontal  velocity  fluctuations  and  the 
near-surface  skewness,  both  in  the  convective  boundary  layer.  In  these  two 
cases  LES  results  are  in  conflict  with  available  experimental  data.  There¬ 
fore,  an  experimental  verification  of  the  numerical  LES  results  is  mandatory 
before  LES  can  be  fully  trusted  as  substitute  for  experiments. 
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Mixing  by  Nonvertical  Shear  in 
the  Stably  Stratified  Ocean 


Sutanu  Sarkar  and  Frank  G.  Jacobitz 

University  of  California  at  San  Diego,  La  Jolla,  CA  92093 

1  Introduction 

Mixing  of  water  masses  of  different  densities,  natural  or  man-made  discharges 
into  the  ocean,  life-supporting  nutrients,  as  well  as  temperature  and  density 
contrasts  is  an  important  determinant  of  local  and  global  features  of  the 
ocean  climate  and  ecosystem.  The  additional  effect  of  stratification  induced 
by  temperature  and  salinity  gradients,  as  well  as  that  of  earth’s  rotation 
must  be  considered  when  mixing  is  studied  in  the  context  of  geophysical  and 
environmental  flows.  At  scales  of  hundred  meters  or  less,  the  effects  of  earth’s 
rotation  may  often  be  neglected.  The  microstructure  at  these  scales  is  thought 
to  be  related  to  a  variety  of  mechanisms  that  include  nonlinearly  breaking 
internal  waves,  shear  layer  instabilities,  convective  instabilities  and  boundary 
layer  turbulence.  Although  diverse  in  their  origin,  mim'ng  processes  in  the 
ocean  often  involve  a  competition  between  shear  and  stable  stratification. 
Therefore,  the  role  of  stably-stratified  shear  flow  in  maintaining  the  ocean 
microstructure  requires  study.  Flows  with  vertical  3hear,  dU / dz ,  have  been 
the  subject  of  numerous  studies;  however,  the  horizontal  shear  component 
dU/dy  has  received  little  to  no  systematic  investigation.  The  influence  of 
horizontal  shear,  dU / dy ,  with  emphasis  on  its  comparison  with  the  oft-studied 
case  of  vertical  shear  is  the  subject  of  this  work. 


Fig.  1.  Sketch  of  (a)  vertically  sheared  flow  and  (b)  horizontally  sheared  flow.  There 
is  a  stable  vertical  density  stratification  in  both  cases 


It  should  be  noted  that,  despite  the  preponderance  of  vertical  shear  in  the 
open  ocean,  there  are  field  measurements,  for  example,  in  coastal  fronts  [1] 
and  straits  [2]  that  indicate  a  significant  dU f  dy  component.  Flows  over  seamounts 
where  greatly  increased  dissipation  rates  have  been  measured  [3]  as  well  as 
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side-  and  bottom-boundary  flows,  most  likely,  have  substantial  horizontal 
shear  components.  To  the  best  of  our  knowledge,  there  have  been  no  theoreti¬ 
cal  or  laboratory  studies  focussed  on  horizontally  sheared,  vertically  stratified 
flow.  However,  laboratory  experiments  relevant  by  virtue  of  the  presence  of 
horizontal  shear  include  the  study  of  a  front  in  a  rotating,  stratified  fluid  [4], 
the  stratified  wake  of  a  sphere  [5],  and  the  stratified  jet  [6]. 

Figures  l(a)-(b)  show  schematics  of  the  case  with  uniform  vertical  mean 
shear,  dUi/dx3  and  uniform  horizontal  mean  shear,  dUi/dxi,  respectively. 
In  both  cases,  the  uniform  mean  stratification  is  vertical  and  stable. 


2  Approach 

Since  little  is  known  from  existing  studies  about  the  influence  of  horizontal 
shear  on  mixing  in  a  vertically  stratified  fluid,  direct  numerical  simulation 
(DNS)  is  an  attractive  approach  that  avoids  the  uncertain  influence  of  tur¬ 
bulence  models  on  the  conclusions.  In  DNS,  all  dynamically  important  scales 
of  motion  are  resolved  which,  given  current  computer  hardware  capabilities, 
limits  the  simulations  to  moderate  Reynolds  number.  In  the  case  of  uniformly 
distorted  flow,  the  appropriate  Reynolds  number  Re\  =  g\/v  is  based  on  the 
rms  velocity,  q  =  y/2K,  and  the  Taylor  microscale,  A  =  ^/5 vq2/e,  and  is 
limited  to  Re\  <  150  in  parametric  investigations.  Here  K  and  e  denote  the 
turbulent  kinetic  energy  and  dissipation  rate,  respectively. 

The  unsteady,  three-dimensional,  Navier-Stokes  equations  with  the  Boussi- 
nesq  approximation  are  numerically  solved.  A  spectral  collocation  method 
is  used  for  the  spatial  discretization  and  a  third-order,  low-storage,  Runge- 
Kutta  scheme  for  the  temporal  advancement  in  order  to  obtain  high  accuracy. 
The  equations  are  solved  in  a  reference  frame  that  moves  with  the  mean  ve¬ 
locity^]  that,  while  introducing  additional  time-dependent  terms  in  the  gov¬ 
erning  equations,  allows  the  use  of  periodic  boundary  conditions  and  thereby 
Fourier  basis  functions.  Up  to  288  x  144  x  144  grid  points  are  used. 

A  series  of  simulations  with  vertical  shear  and  a  second  series  with  hor¬ 
izontal  shear  were  performed  with  the  same  magnitude  of  shear  S  but  with 
different  values  of  the  Brunt-Vaisala  frequency  N  =  (—gSp/po)1^2  obtained 
by  varying  the  mean  stratification  Sp.  Thus,  the  generalized  Richardson  num¬ 
ber  Ri  =  N^/S2  was  varied  in  the  range  0  <  Ri  <  3. 

Although,  the  focus  is  on  contrasting  the  limiting  cases  of  either  horizontal 
or  vertical  shear,  the  effect  of  varying  the  relative  magnitude  of  horizontal 
mean  shear  S2  —  S  cos  9  to  the  vertical  mean  shear  S3  —  S  sin  9  is  also 
ascertained.  The  relative  magnitude  was  varied  by  conducting  a  third  series 
of  simulations  where  the  shear  inclination  angle  9  is  varied  between  9  =  0 
(vertical  shear)  and  9  =  tt/2  (horizontal  shear)  while  keeping  the  magnitude 

of  the  shear  S  -  \J(dUi/dx3)2  +  {<SJ i/dx2)2  constant. 
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3  Results 

The  two  shear  components,  dU [dy  and  dU/dz,  result  in  two  additional  nondi- 
mensional  parameters.  If  two-dimensional  linear  instabilities  in  the  (x,  z) 
plane  are  considered,  the  vertical  Richardson  number  Riv  =  N2/(dU/dy)2  is 
the  only  relevant  nondimensional  parameter  which  would  imply  that  strat¬ 
ification  has  no  effect  on  the  dynamics  when  the  shear  is  horizontal.  DNS 
shows  that  stratification  has  a  substantial  stabilizing  influence  even  when 
the  shear  is  horizontal;  it  is  dear  that  three-dimensional  fluctuations  must 
be  considered. 

The  energetics  of  the  fluctuations  can  be  studied  through  the  transport 
equation  for  turbulent  kinetic  energy,  K,  wherein  the  ratio  of  the  buoyancy 
flux  B  to  the  production  term  P  in  the  turbulence  can  be  written  as 

P  =  fe)  (sin*  S  +  ~vhjvt cos2  0)  (^)  (1) 

after  introducing  the  vertical  and  horizontal  momentum  transport  coefficients 
vv  and  i/ft,  respectively,  as  well  as  the  mass  transport  coefficient,  at-  Equa¬ 
tion  (1)  is  exact  and  implies  that,  two  nondimensional  parameters:  first,  the 
generalized  Richardson  number  Ri  =  IV2 /S2  based  on  the  magnitude  of  the 

shear  5  =  yj {dU i/dx3)2  +  ( dUi/dx-t )2,  and,  second,  the  shear  inclination 
angle  9  are  required. 

In  the  limiting  case  of  purely  horizontal  shear,  dU /dy,  the  relevant  stratification- 
related  parameter  becomes  Ri  =  N2  /  {dU /  dy)2 . 


(a)  a  (b)  s. 


Fig.  2.  Evolution  of  turbulent  kinetic  energy  K  for  various  values  of  Richardson 
number  Ri  in  the  (a)  vertically  sheared  case,  and  (b)  horizontally  sheared  case.  St  is 
time  normalized  with  mean  shear  S.  Dashed  lines  are  the  exponential  approximation 
to  the  solution. 


The  dependence  of  the  turbulence  evolution  on  the  generalized  Richardson 
number,  Ri  =  N2/^,  is  now  discussed.  Figure  2  shows  that,  for  the  same 
value  of  Ri,  the  the  turbulence  is  much  more  energetic  in  the  case  of  horizontal 
shear.  The  critical  Richardson  number,  Ricr,  is  the  value  of  Ri  below  which 
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turbulence  shows  asymptotic  growth  and  above  which  turbulence  decays.  The 
value  of  Ricr  =  1.0  in  the  horizontally  sheared  case  is  an  order  of  magnitude 
larger  than  Ri „.  =  0.15  in  the  vertically  sheared  case.  The  vertical  mass 
transport  is  also  an  order  of  magnitude  larger  when  the  shear  is  horizontal. 
An  explanation  for  these  observations  will  be  given  in  the  full  paper.  In 
addition,  results  on  the  turbulence  mixing  coefficients,  large-  and  small-scale 
anisotropy  and  visualizations  of  mixing  will  be  given  in  the  full  paper. 
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Abstract 

Non-linear  Galerkin  Methods  which  have  been  developed  in  the  context  of  the  long¬ 
time  integration  of  dissipative  evolution  equations  have  been  applied  to  transition 
regimes  of  the  plane  channel  flow  and  the  temporally  growing  free  shear  layer 
problems.  Those  methods  have  been  proposed  in  the  literature  to  approximate  the 
‘inertial  manifolds’  which  are  finite-dimensional  smooth  manifolds  containing  the 
global  attractor,  a  finite-dimensional  manifold  of  complex  structure  and  attracting  all 
the  solutions  at  an  exponential  rate  in  the  phase  space.  These .  manifolds  yield  an 
interaction  law  between  the  small  and  large  scale  components  of  the  flow  and  reduce 
the  number  of  modes  needed  to  describe  the  flow  dynamics.  This  approach  is  thought 
to  be  promising  in  the  context  of  direct  numerical  simulations  of  Navier-Stokes 
equations.  That  work  aims  to  test  the  convergence  and  efficiency  of  different 
Nonlinear  Galerkin  schemes  with  respect  to  each  other  and  with  respect  to  the 
classical  Galerkin  spectral  method  in  the  above  flow  cases. 


Direct  Numerical  Simulation  of  Intrusion  Fronts 


Carlos  Hartel 


Swiss  Federal  Institute  of  Technology 
Institute  of  Fluid  Dynamics 
CH-8092  Zurich,  Switzerland 


Intrusion  fronts  of  heavy  fluid  which  propagate  in  an  environment  of  lighter  fluid  are 
commonly  encountered  in  numerous  geophysical  applications  (see  Simpson,  1997),  well- 
known  examples  being  moving  atmospheric  cold  fronts,  thunderstorm  outflows,  powder- 
snow  avalanches  or  muddy  underflows  in  lakes  or  oceans.  Their  study  is  also  of  interest 
in  the  engineering  sciences  because  of  the  important  role  they  play  in  many  problems 
related  to  industrial  safety  and  environmental  protection  (see  Fannelop,  1994).  A  typical 
scenario  where  intrusion  fronts  may  be  expected  to  form  is  the  accidental  release  of  dense 
(and  possibly  hazardous)  industrial  gases  which  gravity  may  spread  over  relatively  large 
distances  at  an  appreciable  speed.  Extensive  theoretical  and  experimental  studies  were 
devoted  to  intrusion  fronts  in  the  past,  but  very  few  accurate  numerical  simulation 
studies  have  been  presented  so  far.  Most  of  the  numerical  studies  conducted  in  this 
field  have  employed  approximate  methods,  based  e.g.  on  the  shallow-water  equations, 
and  often  utilize  empirical  models  to  account  for  turbulent  transport.  Today,  direct 
numerical  simulations  (DNS)  allow  to  study  the  physics  of  intrusion  fronts  with  much 
greater  accuracy.  DNS  provides  the  full  information  on  the  three-dimensional  and  time- 
dependent  flow  evolution,  and  it  proves  to  be  a  very  powerful  tool  to  address  some  of 
the  intricate  problems  related  to  the  topology  and  stability  characteristics  of  intrusion 
fronts. 

The  talk  will  present  results  from  a  research  project  in  which  the  DNS  technique  is 
applied  for  the  first  time  to  study  fundamental  physical  properties  of  intrusion  fronts. 
Several  prototype  flows  have  been  considered,  one  of  these  being  the  mutual  intrusion 
of  two  fluids  of  different  density  in  a  plane  channel.  Initially  the  two  fluids  are  separated 
by  a  vertical  membrane.  After  the  membrane  has  been  removed  a  heavy-fluid  front  and 
a  light-fluid  front  develop  and  propagate  along  the  lower  and  the  upper  channel  wall, 
respectively.  This  type  of  flow  has  often  been  studied  in  experiments  (see  Grobelbauer 
et  ai,  1993)  and  is  usually  referred  to  as  the  “lock-exchange  problem”.  Distinct  features 
of  lock-exchange  flows  are  the  steep  velocity  and  density  gradients  across  the  heads  of 
the  fronts,  the  thin  boundary  layers  that  form  at  the  walls  and  the  stably  stratified 
interface  between  the  two  fluids  in  the  interior  of  the  channel.  In  Figure  1  the  evolution 
of  a  lock-exchange  flow  at  a  Grashof  number  of  about  10s  (based  on  buoyancy  velocity 
and  channel  half  width)  is  illustrated  by  means  of  isocontours  of  density. 

The  numerical  simulations  are  based  on  the  Boussinesq  equations,  where  density 
variations  are  assumed  to  be  small.  For  the  spatial  discretization  Fourier  expansions 
are  used  in  the  wall-parallel  directions  together  with  a  spectral-element  technique  in 
the  normal  direction.  A  third-order  Runge-Kutta  method  is  employed  for  the  temporal 
discretization  of  the  nonlinear  terms  together  with  a  second-order  accurate  Crank- 
Nicolson  scheme  for  the  viscous  terms  and  the  pressure.  To  validate  the  code  simulations 
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of  Rayleigh-Benard  convection  were  conducted  and  the  results  were  compared  with 
linear  stability  theory  and  reference  data  from  the  literature.  In  ail  cases  an  excellent 
agreement  could  be  obtained.  A  detailed  description  of  the  numerical  scheme  is  given 
in  Hartel  et  al.  (1997). 

Among  other  things,  we  will  discuss  results  of  an  analysis  of  the  pronounced  two- 
dimensional  Kelvin-Helmholtz  type  instability  at  the  interface  between  light  and  heavy 
fluid,  and  also  of  the  three-dimensional  lobe-and-cleft  instability  commonly  observed 
at  the  head  of  propagating  intrusion  fronts.  The  three-dimensional  simulations  had  to 
be  restricted  to  moderate  Grashof  numbers  of  the  order  of  106  due  to  the  excessive 
resolution  requirements.  However,  much  higher  Grashof  numbers  can  be  reached,  if 
strictly  two-dimensional  (2D)  intrusion  fronts  are  considered.  We  have  conducted  2D 
simulations  for  Grashof  numbers  of  up  to  2  •  IQ9  in  order  to  examine  issues  like  the 
dependence  of  the  Froude  number  Ft  of  the  propagating  fronts  (i.e.  the  front  speed 
normalized  by  the  respective  buoyancy  velocity)  on  viscous  forces.  In  Figure  2  the 
relation  between  Ft  and  Gr  is  depicted  for  lock-exchange  simulations  with  no-slip  and 
free-slip  conditions,  respectively,  at  the  top  and  bottom  boundary.  It  is  readily  seen  that 
for  no-slip  walls  the  front  speed  remains  sensitive  to  viscous  effects  over  the  whole  range 
examined.  Moreover  it  can  be  seen  that  our  direct  simulations  are  in  good  agreement 
with  recent  experimental  data. 
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Figure  1:  Numerical  simulation  of  two-dimensional  lock-exchange  flow  at  a  Grashof 
number  of  106.  Shown  are  isocontours  of  density  for  four  successive  time  instants  after 
the  initial  release. 


Figure  2:  Froude  number  Fr  of  the  front  as  a  function  of  Grashof  number  for  lock- 
exchange  flows.  Results  for  no-slip  walls  (solid  line)  and  slip  boundaries.  Symbols 
identify  the  individual  simulations.  The  vertical  bar  gives  the  span  of  results  obtained 
in  recent  lock-exchange  experiments  with  At  and  C02  (Muller  k  Fannelap,  1996). 


Nonlinear  stability  of  the  compressible  attachment-line 
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R.S.  Heeg1  &  B.J.  Geurts,  Faculty  of  Mathematical  Sciences,  University  of  Twente, 
P.O.  Box  217,  7500  AE  Enschede,  The  Netherlands 

The  research  presented  in  this  paper  is  concerned  with  the  stability  properties  of  the 
compressible  attachment-line  boundary  layer.  The  instability  of  flow  near  the  leading 
edge  of  e.g.  a  swept  wing  is  of  great  practical  importance.  This  kind  of  instability  may 
lead  to  growing  disturbances  which  can  be  convected  downstream  and  thus  influence 
the  transition  from  laminar  to  turbulent  flow  around  the  wing.  Thereby  this  process 
may  influence  aerodynamic  properties  to  a  large  extent. 

The  attachment-line  region  consists  of  the  front  part  of  a  wing.  In  figure  I  this  region 
is  shown  schematically  in  order  to  define  the  coordinate  directions.  To  emphasize  the 
far-field  velocity  U  is  at  an  angle  with  respect  to  the  leading  edge  which  is  shown  as 
the  dotted  line.  As  shown  x  denotes  the  chordwise  direction,  y  the  direction  normal  to 


the  wall  and  z  the  spanwise  direction. 

To  address  the  stability  problem  the  flow  is  decomposed  into  a  perturbation  and  a 
basic  flow  part.  The  Navier-Stokes  equations  govern  the  evolution  of  the  perturbation 
around  the  basic  flow.  The  basic  flow  is  the  compressible  counterpart  of  the  swept 
Hiemenz  flow.  The  Navier-Stokes  equations  in  disturbance  form  have  been  used  for 
the  study  the  nonlinear  stability  of  the  attachment-line  boundary  layer. 

As  in  the  incompressible  case  the  linear  eigenmodes  follow  a  sequence  symmetric 
(Si),  antisymmetric  (Al),  symmetric  (S2)  etc.  That  is,  the  most  unstable  mode  is 
symmetric,  the  next  most  unstable  mode  is  antisymmetric,  then  the  next  unstable  mode 
is  symmetric  and  so  on.  In  order  to  study  the  nonlinear  evolution  of  the  disturbances, 
the  linear  perturbations  are  used  as  initial  fields  in  a  number  of  three-dimensional 
temporal  direct  numerical  simulations  (DNS).  These  simulations  employ  high-order 
finite-difference  discretization  in  space  and  implicit  time  integration  using  the  Crank 

1  Present  address:  J.  Huizingalaan  233,  1066  AN  Amsterdam,  The  Netherlands. 


72 


4 


Figure  2:  Growth  rate  as  function  of  time  at  R  =  800,  At  =  0.2,  symmetric  and  general 
simulation. 

Nicolson  scheme.  In  the  z-direction  periodic  boundary  conditions  are  prescribed.  These 
nonlinear  simulations  have  been  validated  using  both  resolution  studies  and  linear 
stability  theory.  The  results  computed  with  the  implicit  time  integration  scheme  have 
also  been  compared  with  results  computed  with  a  conventional  explicit  Runge-Kutta 
scheme  in  order  to  assess  the  time  step  used.  It  was  found  that  the  use  of  the  implicit 
scheme  resulted  in  much  lower  CPU-times  than  the  use  of  the  explicit  scheme  without 
affecting  the  accuracy  of  the  results. 

The  perturbations  could  be  restricted  to  be  symmetric  with  respect  to  the  attach¬ 
ment-line..  This  is  in  fact  an  extension  of  two-dimensional  nonlinear  calculations  of 
previous  researchers  based  upon  the  Gbrtler-Hammerlin  assumption.  A  comparison  has 
been  made  between  the  evolution  of  these  three-dimensional  symmetric  disturbances 
and  the  case  were  no  restriction  on  the  shape  of  the  perturbations  has  been  made.  It 
has  been  found  that,  for  general  disturbances,  the  results  start  to  deviate  significantly 
from  linear  stability  theory  for  disturbance  levels  of  about  1%.  At  these  disturbance 
levels  the  growth  rate  starts  to  increase  to  a  much  higher  value,  see  figure  2.  The 
behaviour  of  the  it-  and  the  T-disturbances  seems  to  be  entirely  responsible  for  the 
deviation  from  linear  theory.  However,  when  only  symmetric  modes  are  allowed  in  the 
simulations,  the  sudden  increase  of  the  growth  rate  is  absent.  Therefore,  the  interaction 
between  symmetric  and  antisymmetric  modes  is  likely  to  be  responsible  for  the  sudden 
increase  of  the  growth  rate. 

In  the  past  no  such  nonlinear  interactions  between  symmetric  and  antisymmetric 
modes  were  found  for  two  reasons.  A  number  of  researchers  used  only  symmetric 
modes  in  their  nonlinear  simulations,  e.g.  Hall  et  al.  (1986)  and  Theofilis  (1998).  Other 
researchers  did  use  a  general  model  for  the  perturbations,  but  used  extremely  small 
disturbance  amplitudes,  e.g.  Spalart  (1988)  and  Joslin  (1995).  Spalart  also  studied 
turbulent  incompressible  attachment-line  flow,  but  the  transition  to  turbulence  was 
not  considered  by  him. 
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DNS/LES  of  Turbulent  Flow  in  a  Square  Duct:  A  Priori 
Evaluation  of  Subgrid  Models 

Peter  L.  O’Sullivan  and  Sedat  Biringen, 

Department  of  Aerospace  Engineering, 

University  of  Colorado, 

Boulder,  CO  80309-0429 

and 

Asmund  Huser, 

Det  Norske  Veritas, 

.  Oslo,  Norway. 

The  direct  simulation  code  for  incompressible  flow  in  a  straight  square  duct  developed  by  Huser 
and  Biringen  [1]  has  been  employed  to  generate  a  high-resolution  (101  x  101  x  128)  database  at 
Rer  =  600  consisting  of  80  total  flow  realizations  which  are  decorrelated  in  time.  The  data  were 
then  post-processed  to  perform  a  priori  testing  of  two  dynamic  SGS  turbulence  models,  namely  the 
dynamic  Smagorinsky  SGS  model  (DSM)  and  the  dynamic  two-parameter  mixed  model  (DTM). 

For  the  duct  flow  there  are  two  inhomogeneous  directions  which  lead  to  “non-canonical”  over¬ 
lapping  turbulent  boundary  layers  presenting  a  challenge  to  more  conventional  approaches  to  large- 
eddy  simulation.  The  duct  can  be  thought  of  as  a  testing  ground  for  developing  subgrid  models 
capable  of  capturing  turbulent  comer  flow  dynamics.  The  two  point  correlation  tensor  also  furnishes 
us  with  estimates  of  the  size  of  large-scale  flow  features,  in  particular  quasi-streamwise  vortices  in 
the  comer  region  of  the  duct.  This  information  serves  as  a  guideline  in  the  minimum  grid-spacing 
required  for  large-eddy  simulations  of  comer  flows  in  which  the  “large  eddies”  are  indeed  captured 
by  the  computer  simulation.  The  mean  diameter  of  streamwise  vortices  appears  to  scale  on  inner 
variables  in  the  near-wall  region  although  this  fact  has  not  been  established  rigorously  for  comer 
flows  of  this  type.  If  inner  scaling  is  indeed  appropriate  then  there  will  be  a  stringent  grid  resolution 
requirement  in  high  Reynolds  number  duct  flows  in  order  to  adequately  capture  the  “smallest  large 
scales”  which  are  dynamically  significant  in  the  comers. 

The  individual  flow  fields  were  also  used  to  perform  a  priori  tests  of  both  the  DSM  and  the  DTM. 
The  success  of  these  models  for  large  eddy  simulations  in  comer  flows  (with  secondary  flows  of  the 
second  kind)  is  as  yet  unexplored.  The  mixed  models  studied  recently  by  Zang,  Street  and  Koseff  [2], 
Liu,  Meneveau  and  Katz  [3]  and  Salvetti  and  Banerjee  [4]  have  shown  better  correlation  with  the 
true  subgrid  stresses.  The  mixed  model  reduces  the  magnitude  of  the  Smagorinsky  parameter  and 
also  generally  makes  it  positive  (enhanced  effective  viscosity)  and  hence  alleviates  the  drawbacks 
(viz.  mathematical  inconsistency,  negative  total  viscosity  and  numerical  instability)  of  the  DSM. 
The  main  reason  for  this  is  the  applicability  of  Bardina’s  scale  similarity  model  to  the  cross  terms 
in  the  subgrid  scale  stresses.  The  scale-similar  part  of  the  model  also  permits  non-alignment  of  the 
principal  axes  of  the  SGS  stresses  with  those  of  the  resolved  strain  rate  tensor  which  is  not  the  case 
in  any  of  the  Smagorinsky-type  models. 

Our  evaluations  currently  consist  of  filtering  individual  DNS  realizations  in  the  ensemble  in 
just  the  streamwise  direction  of  the  flow  which  is  treated  via  periodic  boundary  conditions  in  the 
simulation  code.  Filtering  in  the  inhomogeneous  directions  is  not  well  defined  since  filtering  and 
spatial  differentiation  do  not  commute,  but  rather  introduce  second  order  differencing  errors  whose 
role  is  difficult  to  assess.  Hence,  we  retain  the  high  spatial  resolution  in  the  wall-normal  directions. 
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(The  extension,  to  non-uniform  grid  filters  and  application  to  wall  bounded  spatial  filtering  is  under 
investigation  currently.)  We  base  our  model  test3  on  individual  realizations,  the  SGS  dissipation 
term,  ensemble-average  behavior  and  also  on  the  correlation  coefficient  of  modelled  to  exact  SGS 
dissipation,  p(e).  For  the  DSM  we  find  that  p{e)  «  0.1  —  0.2  at  best  in  agreement  with  Clark, 
Ferziger  and  Reynolds  [5].  The  correlation  in  the  comer  region  is  largely  similar  to  that  in  the  wall 
bisector  region.  For  the  DTM  we  find  that  p(e)  w  0.8  —  0.9  indicating  the  (expected)  improvement 
gained  in  using  a  mixed  model.  The  correlation  is  not  diminished  in  the  comer  regions  of  the  flow. 
The  accompanying  figure  illustrates  this  result. 


Figure  1:  Correlation  coefficient  for  SGS  dissipation  along  the  comer  bisector  for  the  DTM. 

A  favorable  outcome  of  a  priori  model  testing/evaluation  does  not  necessarily  provide  a  suffi¬ 
cient  condition  for  a  successful  LES  using  a  given  model  but  it  does  provide  a  necessary  condition. 
A  model  which  performs  poorly  in  a  priori  tests  stands  little  chance  of  providing  good  temporal 
fidelity  in  LES  even  though  it  may  provide  some  measure  of  success  in  the  mean  (the  success  of 
early  LES  using  the  Smagorinsky  model  are  a  case  in  point).  The  ultimate  test  of  a  given  model  is 
in  its  success  when  utilized  in  a  LES  of  high  Reynolds  number  turbulence. 

We  will  present  results  of  this  work  and  discuss  prospects  for  successful  large  eddy  simulations 
in  similar  complex  flows  of  this  kind. 
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SHALLOW  WATER  MODEL  THE  BOSPHORUS  CURRENT 

HAL  UK  ORS 
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INTRODUCTION 

Bosphorus  strait  is  a  quite  narrow  channel  with  a  mean  width  of  about  a  kilometer  and  length 
of  about  30  km  connecting  the  Black  Sea  to  the  internal  sea  of  Marmara  which  is  considered 
part  of  the  Mediterranean.  The  maximum  depth  is  about  110  meter.  Since  tidal  effects  are 
almost  non-existent  on  the  Bosphorus,  measurements  show  that  tidal  elevation  is  under  4  cm, 
the  flow  is  driven  mainly  by  the  pressure  gradient  in  the  Black  Sea-  Marmara  Sea  direction. 
This  gives  rise  to  a  current  with  a  speed  of  about  4  km/hour  in  the  upper  levels.  The  gradient 
of  the  density  due  to  the  salinity  differences  between  the  Mediterranean  water  of  the  Marmara 
Sea  and  the  less  saline  Black  Sea  water  gives  rise  to  a  countercurrent  in  the  lower  levels  of 
the  Bosphorus.  The  complex  geometry  gives  also  rise  to  local  currents  of  different  directions 
and  the  true  representation  of  the  flow  -field  would  certainly  necessitate  a  three  dimensional 
model  with  a  very  fine  mesh.  However,  surface  currents  are  dominant  and  a  2-D  shallow 
water  model  necessitating  no  more  than  a  desktop  computer  can  be  expected  to  represent  the 
flow  field  with  reasonable  accuracy. 

Application  of  the  Galerkin  finite  element  method  to  the  solution  of  the  shallow  water  equations 
could  be  dated  back  to  the  works  of  Taylor  and  Davis  0.  Norton,  King  and  Iceman  ()  and 
Connor  and  WangO-  These  early  researchers  applied  standard  Galerkin  method  to  discretize 
the  shallow  water  equations  in  space  and  used  Runge-Kutta  or  predictor-corrector  methods  to 
integrate  in  time  the  resulting  system  of  ordinary  differential  equations.  Considering  then 
available  computer  storage  capacities  they  used  highly  coarse  meshes  with  element  sizes  of  the 
order  of  kilometers.  The  well  known  stability  problems  for  standard  Galerkin  methods  were 
circumvented  by  taking  extremely  large  eddy  viscosity  coefficients  or  tuning  the  bottom  friction 
coefficient  for  damping.  For  example  in  the  work  of  Norton  etal  the  diffusivity  coefficient  was 
assumed  to  take  the  value  of  10000  m2/s,  whereas  the  accepted  value  is  around  10  m2/s.  One 
had  to  wait  for  the  development  of  upwinding  techniques  in  the  context  of  finite  elements  for 
the  solution  of  the  convection-difEusion  equation  without  excessive  dissipation 

Shortcomings  of  the  shallow  water  modelling  are  that  the  velocities  computed  represent  simply 
the  height  averaged  values  and  gradients  in  the  vertical  direction  are  lost.  In  cases  where  these 
gradients  are  large  shallow  water  solutions  are  very  difficult  to  compare  with  measured  values. 
Also  it  is  commonly  known  that  surface  waves  in  coastal  seas  generate  a  considerable  high 
velocity  current  in  a  direction  toward  or  outward  the  coast.  This  difficulty  is  usually 
circumvented  by  considering  that  the  depth  is  always  larger  than  a  minimum  value,  which  in 
our  case  is  assumed  to  be  5  meters.  This  assumption  is  quite  correct  in  the  case  of  Bosphorus 
since  both  coasts  can  be  considered  to  form  a  vertically  walled  channel.  In  case  there  is 
beach  special  precautions  need  to  be  taken  as  is  shown  in  (Kawahara,  1984) 

Shallow  water  equations  are  a  good  approximation  to  Navier-Stokes  equations  in  cases  the 
vertical  acceleration  of  the  fluid  can  be  neglected.  This  helps  in  circumventing  the  difficulty 
associated  with  the  computation  of  the  pressure  term  in  the  Navier-Stokes  equation  since  the 
pressure  is  now  hydrostatic  and  easily  computed.  Furthermore  the  problem  is  no  more  three 
dimensional  and  the  number  of  operations  and  storage  requirements  are  at  least  an  order  of 
magnitude  smaller.  Since  the  bathymetry  is  also  taken  into  account,  shallow  water  modeling  is., 
actually  more  than  a  2-dimensional  model  of  the  flow  field  and  can  be  classified  as  a  quasi 


77 


three  dimensional  model.  Finite  element  method  is  the  natural  choice  in  the  derivation  of  the 
discrete  system  of  equations  modeling  phenomena  in  complex  geometries  allowing  us  to  use  an 
unstructured  mesh.  Inclusion  of  the  viscous  effects  by  the  use  of  lateral  eddy  viscosity  makes 
the  system  of  equations  we  use  different  from  the  classical  shallow  water  equations  where  no 
viscosity  term  exists  and  the  type  of  the  equation  is  quasi-hyperbolic.  The  type  of  the  equation 
system  we  will  be  dealing  with  is  incompletely  parabolic  (Gustafsson  and  Sundstrom)  . 

It  is  well-known  that  the  structure  of  the  linear  system  of  equations  has  long  been  a  dominant 
issue  in  the  choice  of  solver  routine  and  usually  effort  is  made  that  the  bandwidth  of  the  final 
matrix  is  as  small  as  possible  to  minimize  the  storage  requirement.  However  the  development 
of  the  matrix-free  solver  routines  based  on  the  Generalized  Minimal  Residual  (GMRES) 
method  in  recent  years  (Saad)  eliminated  this  difficulty.  The  essence  of  the  solver  routine 
consists  in  writing  simply  element  level  matrices  and  collecting  them  in  a  three  indices  array 
without  assembly.  Application  of  the  essential  boundary  conditions  at  the  element  level  further 
reduces  the  number  of  the  equations  to  be  solved.  This  information  together  with  the 
connectivity  table  is  sent  to  the  iterative  solver.  The  advantage  of  this  method  is  to  reduce  the 
memory  requirements  considerably  and  therefore  give  the  possibility  of  attacking  large  scale 
problems  on  a  desktop  computer  .  The  code  is  written  with  consideration  of  a  parallel 
implementation  but  the  results  presented  herein  were  obtained  on  a  single  processor 
workstation. 

As  to  the  time-integration,  trapezoidal  rule  of  time  integration  is  employed.  Trapezoidal  rule 
includes  both  implicit  and  explicit  time  integration  schemes  but  for  the  flexibility  of  the  choice 
of  the  timestep  without  violating  the  CFL  condition  implicit  scheme  is  preferred.  The  use  of  the 
explicit  scheme  is  certainly  more  efficient  computationally  but  the  drawback  is  that  CFL 
condition  restricts  the  size  of  the  timestep  unless  a  specially  designed  mesh  is  used  where  it  is 
made  sure  that,  based  on  the  surface  gravity  wave  speed  -Jgh  for  a  preassigned  time  step, 
CFL  condition  is  not  violated.  (Kashiyama) 

GOVERNING  EQUATIONS 

Averaging  the  velocity  field  in  the  vertical  direction  along  which  fluid  acceleration  is 
neglected  and  therefore  pressure  is  assumed  to  be  simply  hydrostatic,  Navier-Stokes  equations 
are  transformed  into  the  following  set  of  equations, 

~  +  {(/z  +  4')w,.},f  =  0,  (1) 

^T  +  uJuiJ  =  -*(“>•/  +  "yv)  +  r<  -  +/  (2) 

cC 

—  ^C.-kC^s  (3) 

where  represents  the  free  surface  elevation  from  the  mean  water  level,  h  local  mean  water 
depth,  u{  vertically  averaged  horizontal  velocity  components,  C  the  contaminant 
concentration,  e  eddy  viscosity  coefficient,  g  acceleration  due  to  gravity,  r„  wind  stress 
component  at  the  free  surface,  rw  bottom  friction  stress  component,  f  Coriolis  force 
component  and  s  the  contaminant  source  term.  With  this  formulation  the  four  dependent 
variables  of  the  initial  problem,  i.e  ut  for  /=  1,2,3,  p  are  replaced  by  the  three  dependent 
variables  a,  for  /=  1,2,  and  £\  The  difficulty  concerning  the  time  dependent  nature  of  the  free 
surface  is  taken  care  of  by  computing  its  elevation  at  every  timestep  and  using  it  in  computing 
the  local  pressure  gradient  in  conjunction  with  the  local  mean  water  level  h.  The  surface  stress 
term  due  to  the  wind  is  computed  as, 
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(4) 


%  =rPairW,ylW^k 

where  via  is  the  wind  velocity  component,  air  density  and  y  a  friction  coefficient  with  a 
value  of  0.0026  in  SI  system  of  units.  The  bottom  friction  term  is  assumed  to  be  given  by  the 
formula, 

/>  nzg  , - 

-~tfFTuAu*u'<  (5) 

where  n  is  the  Manning  coefficient  with  a  value  of  0.030  in  the  SI  system.  The  contaminant 
diffusion  coefficient  K  is  assumed  to  take  the  value  of  5  m2/s. 

As  to  the  mathematical  classification  of  these  equations,  the  continuity  equation  is  a  first  order 
hyperbolic  equation,  whereas  the  momentum  and  transport  equations  are  second  order  , 
incompletely  parabolic. (Gustaffson) .  Boundary  conditions  appropriate  for  the  well-posedness 
of  these  equations  have  been  subject  to  a  vast  literature  and  according  to  Abraham  et.  al.() 
they  are  as  follows: 

For  subcritical  flow,  which  corresponds  to  our  case,  hyperbolic  equation  requires  2  inflow  and 
1  outflow  boundary  condition,  whereas  the  incompletely  parabolic  equations  require  2 
boundary  conditions  on  land  boundaries,  3  inflow  and  2  outflow  boundary  conditions  at  open 
(sea)  boundaries. 

We  will  see  below  how  these  conditions  translate  into  the  weak  formulation. 

The  simply  connected  open  domain  C2  c  R~  over  which  the  governing  equations  are  valid  is 
closed  by  the  boundary  r  —  rx^jP2  where  Fx  denotes  the  land  type  boundary  and  /',  ocean 
type  boundary,  with  i"j  r\Tz  =  0.  On  land  type  boundary  there  is  no  effluent  flux  and  the 
boundary  condition  is  reduced  to  the  adherence  condition  for  viscous  fluid.  However  since  we 
are  not  interested  in  resolving  the  boundary  layer  and  our  mesh  size  is  extremely  large 
compared  to  the  boundary  layer  thickness  we  will  relax  this  condition  by  cancelling  the  normal 
component  of  the  velocity.(Gray) 

uini  =  0  on  r2  (6) 

On  ocean  type  boundary ,  weak  formulation  gives  rise  to  a  boundary  condition  of  the  form 

+  S(U‘J  +  uij)nj  =  l  on  rx  (7) 

Since  usually  we  have  no  way  of  measuring  the  shear  stress  term  along  the  open  boundary  we 
will  enforce  this  condition  after  neglecting  the  viscosity  term  and  setting  the  tangential 
component  of  the  velocity  at  the  inlet  equal  to  zero.  Therefore  on  the  ocean  boundary  ,  the 
pressure  being  simply  hydrostatic  the  following  condition  will  be  enforced, 

C  =  £on  rx.  (8) 

uinj  =  ~u  on  /j  for  the  inlet  part  only. 

In  the  following  we  prefer  to  work  with  the  total  height  H  =  h+  £ . 

FINITE  ELEMENT  DISCRETIZATION 

Flow  field  being  divided  into  triangular  elements,  the  dependent  variables  are  represented  over 
each  element  by  the  linear  interpolating  functions  0, .  This  is  equivalent  into  discretizing  the 
dependent  variables  in  terms  of  their  nodal  values  as. 
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(9) 


H  =  Ha®a .  ui  =  «,<A ,  and  C  =  cad>a 

where  repeated  indices  mean  summation  unless  otherwise  stated. 

Considering  the  size  of  the  elements  used,  side  length  between  about  100  m  and  300  m,  and 
the  thickness  of  the  boundary  layer  on  the  coast  it  might  have  been  more  accurate  to  use  no 
flux  boundary  condition  along  the  coast  but  with  the  primitive  variables  formulation  of  the 
governing  equations  this  type  of  boundary  condition  is  difficult  to  be  implemented  on  a 
boundary  not  necessarily  aligned  with  the  main  axis.  Stream  function-vorticity  formulation  is 
certainly  more  suitable  in  this  sense  but  -then  linear  triangular  elements  are  unable  to  take 
account  of  the  higher  order  derivatives  originating  from  the  viscous  terms.  Higher  order 
elements  would  be  needed  and  for  the  same  mesh  size  that  would  at  least  double  the  number 
equations  to  be  solved  in  each  time  step.  The  approach  mentioned  by  Gray  ()f°r  the 
implementation  of  slip  boundary  conditions  by  rotation  of  the  element  matrices  corresponding 
to  boundary  elements  is  a  method  operational  within  the  context  of  linear  elements  but  has  not 
been  attempted  in  this  work. 

Initially  it  is  assumed  that  the  deviation  from  the  mean  water  level  is  zero  everywhere  except  on 
the  boundary  at  the  Black  Sea  where  a  constant  deviation  of  0.30  m,  based  on  measurements 
between  Black  Sea  and  Marmara  Sea  levels,  is  assumed  and  maintained  during  the 
computation.  This  difference  in  water  level  creates  a  gravity  travelling  along  the  strait  and 
eventually  sets  up  the  velocity  field.  The  time  the  gravity  wave  takes  to  traverse  the  strait  is 
about  3000  s  which  corresponds  to  a  wave  speed  of  10  m/s.  This  is  the  wave  speed 
corresponding  to  a  constant  depth  of  approximately  10  m..  The  velocity  field  of  interest 
corresponding  to  a  quasi  steady  state  is  obtained  after  t=3000  s.  This  velocity  distribution 
could  afterwards  be  used  as  initial  velocity  distribution  computations  with  wind  loadings  and 
the  contaminant  transport  equation  solution.  Most  common  winds  in  the  Bosphorus  are 
northerly  with  eventual  southerly  winds  in  winters.  In  the  computations  presented  herein  the 
wind  speed  is  taken  to  be  about  20  km/h  and  the  direction  is  northerly. 
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